{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An old FiPy solution to 1D BL\n",
    "I'm not really sure if the result is correct.\n",
    "The results is indeed incorrect. See the updated code a few cells below that might work slightly better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAEICAYAAAC5yopxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZo0lEQVR4nO3deXRc5Z3m8e9TJVmyjc1i2TR4iUxjCA5eIGLLQhsYCFsH+iSZZnUgdHyAMEn6NGlg+gw45AxNn2aabg4EjydDSDNhmQCHOMYNCQFPIJAOcsJmwGCwwcIJyKYx3m1Jv/mjSkUhZKusuqLqXj+fc4R0l6r7U6F6/N73vnVfRQRmZgC5WhdgZvXDgWBmJQ4EMytxIJhZiQPBzEocCGZW4kCwj42kcyX9vNZ12I7J4xB2T5ICmBIRy4fo+VuBFUBjRHQNxTEseW4h2KBIyte6BkueAyEDJF0h6S1J6yUtk3SCpCMlPSXpPUl/kHSzpGHF/X9VfOizkjZI+ktJF0h6os/zhqQDiz/fLulWSYskbQSOk3SapN9Lel/SKklzyx7ee4z3isc4pu8xJH1G0tOS1hW/f6Zs22JJ35P06+Lv9XNJLUPw8lkZB0LKSToYuAw4IiJGAV8AVgLdwF8DLcAxwAnApQARcWzx4TMiYo+IuKfCw50D/HdgFPAEsBGYDewFnAZcIunM4r69x9ireIyn+tS9D/AgcBMwBvgn4EFJY/oc70JgHDAMuLzCOm2QHAjp1w00AVMlNUbEyoh4LSKWRMRvIqIrIlYC/xP4syqP9dOI+HVE9ETElohYHBHPF5efA+7ahWOcBrwaEXcUa7wLeBn487J9fhgRr0TEZuD/AjOrrN8G4EBIuWKn4LeBucA7ku6WtL+kgyQtlPRHSe8D11FoLVRjVfmCpKMkPSapU9I64OJdOMb+wBt91r0BjC9b/mPZz5uAPXaxXttFDoQMiIg7I+JzwCeAAP4BuJXCv7hTImI08F8B7eRpNgIjehck/Ul/h+qzfCewAJgYEXsC88qOMdDlq9XFestNAt4a4HE2hBwIKSfpYEnHS2oCtgCbKZxGjALeBzZI+iRwSZ+Hvg0cULb8LPApSTMlNVNocQxkFPBuRGyRdCSFc/5enUBPn2OUWwQcJOkcSQ2S/hKYCiys4Lg2RBwI6dcEXA+sodDEHkehNXA5hTfoeuB/AX07DucCPypehfjPEfEKcC3wCPAqhU7DgVwKXCtpPXA1hfN8ACJiE4UOyF8Xj3F0+QMjYi1wOvA3wFrgb4HTI2JN5b+6Jc0Dk8ysxC0EMytxIJhZiQPBzEocCGZW0lCrA7e0tERra2utDm+221qyZMmaiBjb37aaBUJrayvt7e21OrzZbktS3xGiJT5lMLMSB4KZlTgQzKykZn0I/dm+fTsdHR1s2bKl1qVkSnNzMxMmTKCxsbHWpVidq6tA6OjoYNSoUbS2tiLt7IN5VqmIYO3atXR0dDB58uRal2N1rq5OGbZs2cKYMWMcBgmSxJgxY9zqsorUVSAADoMh4NfUKlV3gWBmteNAqMLtt9/O6tWrE3u+lStXcuedd5aW29vb+eY3v5nY85sNxIFQhcEEQlfXjucs6RsIbW1t3HTTTYOuz2xXORD62LhxI6eddhozZszg0EMP5Z577uHaa6/liCOO4NBDD2XOnDlEBPfeey/t7e2ce+65zJw5k82bN9Pa2sqaNYUb/rS3tzNr1iwA5s6dy5w5czjppJOYPXs2K1eu5POf/zyHH344hx9+OE8++SQAV155JY8//jgzZ87kxhtvZPHixZx++ukAvPvuu5x55plMnz6do48+mueee6703F/72teYNWsWBxxwgAPEqlJXlx3LffdnS3lx9fuJPufU/UdzzZ9/aqf7PPTQQ+y///48+OCDAKxbt44TTzyRq6++GoDzzz+fhQsX8uUvf5mbb76ZG264gba2tgGPvWTJEp544gmGDx/Opk2b+MUvfkFzczOvvvoqZ599Nu3t7Vx//fXccMMNLFxYuK3g4sWLS4+/5pprOOyww3jggQd49NFHmT17Ns888wwAL7/8Mo899hjr16/n4IMP5pJLLvGYAxsUtxD6mDZtGo888ghXXHEFjz/+OHvuuSePPfYYRx11FNOmTePRRx9l6dKlu/y8X/ziFxk+fDhQGID19a9/nWnTpvGVr3yFF198ccDHP/HEE5x//vkAHH/88axdu5Z169YBcNppp9HU1ERLSwvjxo3j7bff3uX6zKCCFoKk2yjcDPOdiDi0n+0C/gU4lcK98y+IiN9VW9hA/5IPlYMOOoglS5awaNEirrrqKk466SRuueUW2tvbmThxInPnzt3hNf2GhgZ6enoAPrLPyJEjSz/feOON7Lvvvjz77LP09PTQ3Nw8YF393fuy93JiU1NTaV0+n99pP4XZzlTSQrgdOHkn208BphS/5lCYDyC1Vq9ezYgRIzjvvPO4/PLL+d3vCtnW0tLChg0buPfee0v7jho1ivXr15eWW1tbWbJkCQD33XffDo+xbt069ttvP3K5HHfccQfd3d39Pl+5Y489lh//+MdA4VSipaWF0aNHV/fLmvUxYAshIn5VnNp7R84A/jUK/4T9RtJekvaLiD9UW9y2rh62dnV/aN1Hh9joQysHGoKj4n9yEg05kc/pQwN3nn/+eb7zne+Qy+VobGzk1ltv5YEHHmDatGm0trZyxBFHlPa94IILuPjiixk+fDhPPfUU11xzDRdddBHXXXcdRx111A5ruPTSS/nSl77ET37yE4477rhS62H69Ok0NDQwY8YMLrjgAg477LDSY+bOncuFF17I9OnTGTFiBD/60Y8G+E3Ndl1Ft2EvBsLCHZwyLASuj4gnisu/BK6IiI/c/UTSHAqtCCZNmvTpN9748H0aXnrpJQ455JDS8toNW3nrvc278OvsOiGaGnOMaMyzR3MDo5sbyeWyN7Kv72truy9JSyKi357wJK4y9Pfu6TdlImI+MB+gra1twCQaPbyR5sb8jp84PnqonT1p+e7dEXR1B109PWzZ3sO6Ldt5d9M2GnI5xo1uYszIYR7ya7udJAKhA5hYtjyBwrx9VWvM52jMfzwXQiKCDVu76Fy/ldXvbWbdpu1MGjPiYzu+WT1I4q99ATBbBUcD66rpP6jVTFKSGNXcyOSWkUzcewSbt3fzWucGtvXpw0gjz85llarksuNdwCygRVIHcA3QCBAR8yhM2nkqsJzCZccLB1tMc3Mza9eurelHoCWx98hhNDXmWLFmIyvWbOJPx46kIaUthd77IVRyadOskqsMZw+wPYBvJFHMhAkT6OjooLOzM4mnq9r2rh7+sGErnatyjNmjaeAH1KneOyaZDaSuhi43NjbW3V19fvjrFXz3vhe57i+mcc5Rk2pdjtmQSmc7+GP01WNaOeaAMfz9v73E2g1ba12O2ZByIAwglxPXnvEpNm3r5sZHXql1OWZDyoFQgSn7juLsIydyz9OrhnyglFktORAqdMmsAwGYt/i1GldiNnQcCBUav9dwzpw5nnuXdPD+lu21LsdsSDgQdsHsY1rZvL2b+5d01LoUsyHhQNgF0ybsyYwJe3LXb1fVuhSzIeFA2EVf+vQElr29nmV/7P++BWZp5kDYRaccuh85wc+eTe7262b1woGwi8aOauIzf9rCouervv+LWd1xIAzCCYeM4/U1G3lj7cZal2KWKAfCIMw6eBwAi5fVx4ewzJLiQBiEyS0j+cSYESxe9k6tSzFLlANhkD53YAu/XfEu3T2++YhlhwNhkNpa92bjtm5ffrRMcSAMUtsn9gFgyZv/UeNKzJLjQBikCXsPZ+yoJn73hgPBssOBMEiS+PSkvfm9WwiWIQ6EKhy07x68+e6mj8wuZZZWDoQqHDB2D3oCVr27qdalmCXCgVCFA8YW5mR8rdMjFi0bHAhVmNxSCITXHQiWEQ6EKoxqbmTsqCZe79xQ61LMEuFAqNIBLSNZscYtBMsGB0KVJuw9wnditsxwIFSpuTHH1q6eWpdhlggHQpWaGvJscyBYRlQUCJJOlrRM0nJJV/azfU9JP5P0rKSlkgY9A3TaDGvIORAsMwYMBEl54BbgFGAqcLakqX12+wbwYkTMoDB1/P+QNCzhWuvSsIYc27p76PHHoC0DKmkhHAksj4jXI2IbcDdwRp99AhglScAewLtAV6KV1qmmhsJLuK3brQRLv0oCYTxQPhFBR3FduZuBQ4DVwPPAtyLiI+8QSXMktUtq7+zMxu3HHAiWJZUEgvpZ17d9/AXgGWB/YCZws6TRH3lQxPyIaIuItrFjx+5ysfVoWG8guB/BMqCSQOgAJpYtT6DQEih3IXB/FCwHVgCfTKbE+jYsX3gJfenRsqCSQHgamCJpcrGj8CxgQZ993gROAJC0L3Aw8HqShdarpka3ECw7GgbaISK6JF0GPAzkgdsiYqmki4vb5wHfA26X9DyFU4wrImLNENZdN4bl84ADwbJhwEAAiIhFwKI+6+aV/bwaOCnZ0tKhtw/BN0mxLPBIxSo1uVPRMsSBUCVfZbAscSBUqXTK4HEIlgEOhCr1njJs3e5AsPRzIFTJIxUtSxwIVfJlR8sSB0KVfNnRssSBUCVfdrQscSBUyZcdLUscCFVyIFiWOBCq1JATOfnTjpYNDoQqSSrdRs0s7RwICRiW941WLRscCAkY1pD3KYNlggMhAU0NOY9DsExwICSgyXMzWEY4EBLgyVosKxwICSicMjgQLP0cCAlwC8GywoGQAI9DsKxwICTA4xAsKxwICWhqyPuyo2WCAyEB7kOwrHAgJMCBYFnhQEhAYz7Htu6+89+apY8DIQH5HPSEA8HSz4GQgLxEd48DwdKvokCQdLKkZZKWS7pyB/vMkvSMpKWS/l+yZda3fC7nQLBMGHCyV0l54BbgRKADeFrSgoh4sWyfvYDvAydHxJuSxg1VwfUon8OBYJlQSQvhSGB5RLweEduAu4Ez+uxzDnB/RLwJEBHvJFtmfcvlRLf7ECwDKgmE8cCqsuWO4rpyBwF7S1osaYmk2f09kaQ5ktoltXd2dg6u4jqUl+hxC8EyoJJAUD/r+v71NwCfBk4DvgD8N0kHfeRBEfMjoi0i2saOHbvLxdarvFsIlhED9iFQaBFMLFueAKzuZ581EbER2CjpV8AM4JVEqqxz+ZyIgJ6eIJfrLz/N0qGSFsLTwBRJkyUNA84CFvTZ56fA5yU1SBoBHAW8lGyp9SuvQgi4lWBpN2ALISK6JF0GPAzkgdsiYqmki4vb50XES5IeAp4DeoAfRMQLQ1l4PeltFXT3BI35GhdjVoVKThmIiEXAoj7r5vVZ/kfgH5MrLT3yxUDwaEVLO49UTEDplMFXGizlHAgJ6D1l6PEHHi3lHAgJaMi5U9GywYGQgN4WQpebCJZyDoQE9PYhOA8s7RwICcgXX0WfMljaORASkCu1EBwIlm4OhATkc77saNngQEhA3lcZLCMcCAlwC8GywoGQAI9UtKxwICQg5xaCZYQDIQGlcQjuQ7CUcyAkwH0IlhUOhAT448+WFQ6EBPQGQpenc7OUcyAkIOdbqFlGOBASkPf9ECwjHAgJ8IebLCscCAnwh5ssKxwICWjIFV5GX3a0tHMgJKCYB3Q5ECzlHAgJ8DgEywoHQgL84SbLCgdCAnJuIVhGOBAS4BaCZYUDIQH+cJNlRUWBIOlkScskLZd05U72O0JSt6QvJ1di/XMgWFYMGAiS8sAtwCnAVOBsSVN3sN8/UJglerfieypaVlTSQjgSWB4Rr0fENuBu4Ix+9vsvwH3AOwnWlwoeqWhZUUkgjAdWlS13FNeVSBoP/AXwoSni+5I0R1K7pPbOzs5drbVu+ZTBsqKSQFA/6/r+5f8zcEVEdO/siSJifkS0RUTb2LFjK62x7pWuMjgPLOUaKtinA5hYtjwBWN1nnzbgbhXeGC3AqZK6IuKBRKqsc/l8bwvBn3+2dKskEJ4GpkiaDLwFnAWcU75DREzu/VnS7cDC3SUMoHwcQo0LMavSgIEQEV2SLqNw9SAP3BYRSyVdXNy+036D3UHvh5s8UtHSrpIWAhGxCFjUZ12/QRARF1RfVrp4pKJlhUcqJsBXGSwrHAgJkITkUwZLPwdCQhpy8g1SLPUcCAnJSR6paKnnQEhIPif3IVjqORASkpf84SZLPQdCQnI5nzJY+jkQEpLPuYVg6edASIj7ECwLHAgJycuBYOnnQEhIoYVQ6yrMquNASEgu55GKln4OhIT4lMGywIGQkJyvMlgGOBAS0pAT3b6HmqWcAyEhOY9UtAxwICQk75GKlgEOhIR4pKJlgQMhITlfZbAMcCAkJJ+TxyFY6jkQEpLPiS5fZbCUcyAkJC+3ECz9HAgJ8acdLQscCAkpjFSsdRVm1XEgJCQvTwdv6edASIhPGSwLHAgJcSBYFlQUCJJOlrRM0nJJV/az/VxJzxW/npQ0I/lS65tHKloWDBgIkvLALcApwFTgbElT++y2AviziJgOfA+Yn3Sh9c4TtVgWVNJCOBJYHhGvR8Q24G7gjPIdIuLJiPiP4uJvgAnJlln/3EKwLKgkEMYDq8qWO4rrduQi4N/62yBpjqR2Se2dnZ2VV5kCvmOSZUElgaB+1vX7ly/pOAqBcEV/2yNifkS0RUTb2LFjK68yBTxRi2VBQwX7dAATy5YnAKv77iRpOvAD4JSIWJtMeenh2Z8tCyppITwNTJE0WdIw4CxgQfkOkiYB9wPnR8QryZdZ/3L+tKNlwIAthIjoknQZ8DCQB26LiKWSLi5unwdcDYwBvi8JoCsi2oau7PrjFoJlQSWnDETEImBRn3Xzyn7+K+Cvki0tXYblc2zr8kwtlm4eqZiQ5sY8Wx0IlnIOhIQ0NeTo7gm6PJ+bpZgDISFNjYWXcotbCZZiDoSENDXkAdi6vbvGlZgNngMhIc3FFoL7ESzNHAgJKbUQHAiWYg6EhDQ1FPsQfMpgKeZASEiTTxksAxwICWl2p6JlgAMhIW4hWBY4EBLS26noPgRLMwdCQno7Fd1CsDRzICSkudGXHS39HAgJ+aCF4FMGSy8HQkI+6ENwC8HSy4GQkA+uMriFYOnlQEhI6ZTBLQRLMQdCQiQxrCHnTkVLNQdCgpoach6HYKnmQEiQb6NmaedASFBTQ86dipZqDoQENbkPwVLOgZCgpoa8P+1oqeZASFBzo1sIlm4OhAQVWggOBEsvB0KCmhrdqWjp5kBIUGEcglsIll4OhAQVxiG4hWDpVVEgSDpZ0jJJyyVd2c92SbqpuP05SYcnX2r982VHS7sBZ3+WlAduAU4EOoCnJS2IiBfLdjsFmFL8Ogq4tfh9t9LUkGf9li7eem8zI4fla12O7Yb2HN6IpEE/vpLp4I8ElkfE6wCS7gbOAMoD4QzgXyMigN9I2kvSfhHxh0FXlkKfPbCFO3/7Jp+9/tFal2K7qRV/f2pVj68kEMYDq8qWO/jov/797TMe+FAgSJoDzAGYNGnSrtZa904+9E94+NvH8sSrnUStizEbhEoCob/2R9+/90r2ISLmA/MB2traMvmeOXDcHhw4bo9al2E2KJV0KnYAE8uWJwCrB7GPmdW5SgLhaWCKpMmShgFnAQv67LMAmF282nA0sG536z8wy4IBTxkiokvSZcDDQB64LSKWSrq4uH0esAg4FVgObAIuHLqSzWyoVNKHQEQsovCmL183r+znAL6RbGlm9nHzSEUzK3EgmFmJA8HMShwIZlbiQDCzEgeCmZU4EMysxIFgZiUOBDMrcSCYWYkDwcxKHAhmVqLC55JqcGCpE3ijgl1bgDVDXE616r3Geq8PXGNSKqnxExExtr8NNQuESklqj4i2WtexM/VeY73XB64xKdXW6FMGMytxIJhZSRoCYX6tC6hAvddY7/WBa0xKVTXWfR+CmX180tBCMLOPiQPBzErqNhAGmmC2FiRNlPSYpJckLZX0reL6fST9QtKrxe9717jOvKTfS1pYj/UVa9pL0r2SXi6+nsfUU52S/rr4//gFSXdJaq51fZJuk/SOpBfK1u2wJklXFd8/yyR9oZJj1GUglE0wewowFThb0tTaVgVAF/A3EXEIcDTwjWJdVwK/jIgpwC+Ly7X0LeClsuV6qw/gX4CHIuKTwAwK9dZFnZLGA98E2iLiUArTD5xVB/XdDpzcZ12/NRX/Ls8CPlV8zPeL76udi4i6+wKOAR4uW74KuKrWdfVT508pzIq9DNivuG4/YFkNa5pQ/MM4HlhYXFc39RVrGA2soNipXba+Lurkg7lK96EwVcFC4KR6qA9oBV4Y6DXr+56hMK/KMQM9f122ENjx5LF1Q1IrcBjw78C+UZypqvh9XO0q45+BvwV6ytbVU30ABwCdwA+LpzY/kDSSOqkzIt4CbgDepDBh8bqI+Hm91NfHjmoa1HuoXgOhoslja0XSHsB9wLcj4v1a19NL0unAOxGxpNa1DKABOBy4NSIOAzZSH6cxABTPw88AJgP7AyMlnVfbqnbZoN5D9RoIdTt5rKRGCmHw44i4v7j6bUn7FbfvB7xTo/I+C3xR0krgbuB4Sf+njurr1QF0RMS/F5fvpRAQ9VLnfwJWRERnRGwH7gc+U0f1ldtRTYN6D9VrIFQywezHTpKA/w28FBH/VLZpAfDV4s9fpdC38LGLiKsiYkJEtFJ4zR6NiPPqpb5eEfFHYJWkg4urTgBepH7qfBM4WtKI4v/zEyh0etZLfeV2VNMC4CxJTZImA1OA3w74bLXotKmw8+RU4BXgNeDval1PsabPUWh2PQc8U/w6FRhDoSPv1eL3feqg1ll80KlYj/XNBNqLr+UDwN71VCfwXeBl4AXgDqCp1vUBd1Ho09hOoQVw0c5qAv6u+P5ZBpxSyTE8dNnMSur1lMHMasCBYGYlDgQzK3EgmFmJA8HMShwIZlbiQDCzkv8PCcNkFBFhFbIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<fipy.viewers.matplotlibViewer.matplotlib1DViewer.Matplotlib1DViewer at 0x7f8bb6751588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAEICAYAAAC5yopxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZo0lEQVR4nO3deXRc5Z3m8e9TJVmyjc1i2TR4iUxjCA5eIGLLQhsYCFsH+iSZZnUgdHyAMEn6NGlg+gw45AxNn2aabg4EjydDSDNhmQCHOMYNCQFPIJAOcsJmwGCwwcIJyKYx3m1Jv/mjSkUhZKusuqLqXj+fc4R0l6r7U6F6/N73vnVfRQRmZgC5WhdgZvXDgWBmJQ4EMytxIJhZiQPBzEocCGZW4kCwj42kcyX9vNZ12I7J4xB2T5ICmBIRy4fo+VuBFUBjRHQNxTEseW4h2KBIyte6BkueAyEDJF0h6S1J6yUtk3SCpCMlPSXpPUl/kHSzpGHF/X9VfOizkjZI+ktJF0h6os/zhqQDiz/fLulWSYskbQSOk3SapN9Lel/SKklzyx7ee4z3isc4pu8xJH1G0tOS1hW/f6Zs22JJ35P06+Lv9XNJLUPw8lkZB0LKSToYuAw4IiJGAV8AVgLdwF8DLcAxwAnApQARcWzx4TMiYo+IuKfCw50D/HdgFPAEsBGYDewFnAZcIunM4r69x9ireIyn+tS9D/AgcBMwBvgn4EFJY/oc70JgHDAMuLzCOm2QHAjp1w00AVMlNUbEyoh4LSKWRMRvIqIrIlYC/xP4syqP9dOI+HVE9ETElohYHBHPF5efA+7ahWOcBrwaEXcUa7wLeBn487J9fhgRr0TEZuD/AjOrrN8G4EBIuWKn4LeBucA7ku6WtL+kgyQtlPRHSe8D11FoLVRjVfmCpKMkPSapU9I64OJdOMb+wBt91r0BjC9b/mPZz5uAPXaxXttFDoQMiIg7I+JzwCeAAP4BuJXCv7hTImI08F8B7eRpNgIjehck/Ul/h+qzfCewAJgYEXsC88qOMdDlq9XFestNAt4a4HE2hBwIKSfpYEnHS2oCtgCbKZxGjALeBzZI+iRwSZ+Hvg0cULb8LPApSTMlNVNocQxkFPBuRGyRdCSFc/5enUBPn2OUWwQcJOkcSQ2S/hKYCiys4Lg2RBwI6dcEXA+sodDEHkehNXA5hTfoeuB/AX07DucCPypehfjPEfEKcC3wCPAqhU7DgVwKXCtpPXA1hfN8ACJiE4UOyF8Xj3F0+QMjYi1wOvA3wFrgb4HTI2JN5b+6Jc0Dk8ysxC0EMytxIJhZiQPBzEocCGZW0lCrA7e0tERra2utDm+221qyZMmaiBjb37aaBUJrayvt7e21OrzZbktS3xGiJT5lMLMSB4KZlTgQzKykZn0I/dm+fTsdHR1s2bKl1qVkSnNzMxMmTKCxsbHWpVidq6tA6OjoYNSoUbS2tiLt7IN5VqmIYO3atXR0dDB58uRal2N1rq5OGbZs2cKYMWMcBgmSxJgxY9zqsorUVSAADoMh4NfUKlV3gWBmteNAqMLtt9/O6tWrE3u+lStXcuedd5aW29vb+eY3v5nY85sNxIFQhcEEQlfXjucs6RsIbW1t3HTTTYOuz2xXORD62LhxI6eddhozZszg0EMP5Z577uHaa6/liCOO4NBDD2XOnDlEBPfeey/t7e2ce+65zJw5k82bN9Pa2sqaNYUb/rS3tzNr1iwA5s6dy5w5czjppJOYPXs2K1eu5POf/zyHH344hx9+OE8++SQAV155JY8//jgzZ87kxhtvZPHixZx++ukAvPvuu5x55plMnz6do48+mueee6703F/72teYNWsWBxxwgAPEqlJXlx3LffdnS3lx9fuJPufU/UdzzZ9/aqf7PPTQQ+y///48+OCDAKxbt44TTzyRq6++GoDzzz+fhQsX8uUvf5mbb76ZG264gba2tgGPvWTJEp544gmGDx/Opk2b+MUvfkFzczOvvvoqZ599Nu3t7Vx//fXccMMNLFxYuK3g4sWLS4+/5pprOOyww3jggQd49NFHmT17Ns888wwAL7/8Mo899hjr16/n4IMP5pJLLvGYAxsUtxD6mDZtGo888ghXXHEFjz/+OHvuuSePPfYYRx11FNOmTePRRx9l6dKlu/y8X/ziFxk+fDhQGID19a9/nWnTpvGVr3yFF198ccDHP/HEE5x//vkAHH/88axdu5Z169YBcNppp9HU1ERLSwvjxo3j7bff3uX6zKCCFoKk2yjcDPOdiDi0n+0C/gU4lcK98y+IiN9VW9hA/5IPlYMOOoglS5awaNEirrrqKk466SRuueUW2tvbmThxInPnzt3hNf2GhgZ6enoAPrLPyJEjSz/feOON7Lvvvjz77LP09PTQ3Nw8YF393fuy93JiU1NTaV0+n99pP4XZzlTSQrgdOHkn208BphS/5lCYDyC1Vq9ezYgRIzjvvPO4/PLL+d3vCtnW0tLChg0buPfee0v7jho1ivXr15eWW1tbWbJkCQD33XffDo+xbt069ttvP3K5HHfccQfd3d39Pl+5Y489lh//+MdA4VSipaWF0aNHV/fLmvUxYAshIn5VnNp7R84A/jUK/4T9RtJekvaLiD9UW9y2rh62dnV/aN1Hh9joQysHGoKj4n9yEg05kc/pQwN3nn/+eb7zne+Qy+VobGzk1ltv5YEHHmDatGm0trZyxBFHlPa94IILuPjiixk+fDhPPfUU11xzDRdddBHXXXcdRx111A5ruPTSS/nSl77ET37yE4477rhS62H69Ok0NDQwY8YMLrjgAg477LDSY+bOncuFF17I9OnTGTFiBD/60Y8G+E3Ndl1Ft2EvBsLCHZwyLASuj4gnisu/BK6IiI/c/UTSHAqtCCZNmvTpN9748H0aXnrpJQ455JDS8toNW3nrvc278OvsOiGaGnOMaMyzR3MDo5sbyeWyN7Kv72truy9JSyKi357wJK4y9Pfu6TdlImI+MB+gra1twCQaPbyR5sb8jp84PnqonT1p+e7dEXR1B109PWzZ3sO6Ldt5d9M2GnI5xo1uYszIYR7ya7udJAKhA5hYtjyBwrx9VWvM52jMfzwXQiKCDVu76Fy/ldXvbWbdpu1MGjPiYzu+WT1I4q99ATBbBUcD66rpP6jVTFKSGNXcyOSWkUzcewSbt3fzWucGtvXpw0gjz85llarksuNdwCygRVIHcA3QCBAR8yhM2nkqsJzCZccLB1tMc3Mza9eurelHoCWx98hhNDXmWLFmIyvWbOJPx46kIaUthd77IVRyadOskqsMZw+wPYBvJFHMhAkT6OjooLOzM4mnq9r2rh7+sGErnatyjNmjaeAH1KneOyaZDaSuhi43NjbW3V19fvjrFXz3vhe57i+mcc5Rk2pdjtmQSmc7+GP01WNaOeaAMfz9v73E2g1ba12O2ZByIAwglxPXnvEpNm3r5sZHXql1OWZDyoFQgSn7juLsIydyz9OrhnyglFktORAqdMmsAwGYt/i1GldiNnQcCBUav9dwzpw5nnuXdPD+lu21LsdsSDgQdsHsY1rZvL2b+5d01LoUsyHhQNgF0ybsyYwJe3LXb1fVuhSzIeFA2EVf+vQElr29nmV/7P++BWZp5kDYRaccuh85wc+eTe7262b1woGwi8aOauIzf9rCouervv+LWd1xIAzCCYeM4/U1G3lj7cZal2KWKAfCIMw6eBwAi5fVx4ewzJLiQBiEyS0j+cSYESxe9k6tSzFLlANhkD53YAu/XfEu3T2++YhlhwNhkNpa92bjtm5ffrRMcSAMUtsn9gFgyZv/UeNKzJLjQBikCXsPZ+yoJn73hgPBssOBMEiS+PSkvfm9WwiWIQ6EKhy07x68+e6mj8wuZZZWDoQqHDB2D3oCVr27qdalmCXCgVCFA8YW5mR8rdMjFi0bHAhVmNxSCITXHQiWEQ6EKoxqbmTsqCZe79xQ61LMEuFAqNIBLSNZscYtBMsGB0KVJuw9wnditsxwIFSpuTHH1q6eWpdhlggHQpWaGvJscyBYRlQUCJJOlrRM0nJJV/azfU9JP5P0rKSlkgY9A3TaDGvIORAsMwYMBEl54BbgFGAqcLakqX12+wbwYkTMoDB1/P+QNCzhWuvSsIYc27p76PHHoC0DKmkhHAksj4jXI2IbcDdwRp99AhglScAewLtAV6KV1qmmhsJLuK3brQRLv0oCYTxQPhFBR3FduZuBQ4DVwPPAtyLiI+8QSXMktUtq7+zMxu3HHAiWJZUEgvpZ17d9/AXgGWB/YCZws6TRH3lQxPyIaIuItrFjx+5ysfVoWG8guB/BMqCSQOgAJpYtT6DQEih3IXB/FCwHVgCfTKbE+jYsX3gJfenRsqCSQHgamCJpcrGj8CxgQZ993gROAJC0L3Aw8HqShdarpka3ECw7GgbaISK6JF0GPAzkgdsiYqmki4vb5wHfA26X9DyFU4wrImLNENZdN4bl84ADwbJhwEAAiIhFwKI+6+aV/bwaOCnZ0tKhtw/BN0mxLPBIxSo1uVPRMsSBUCVfZbAscSBUqXTK4HEIlgEOhCr1njJs3e5AsPRzIFTJIxUtSxwIVfJlR8sSB0KVfNnRssSBUCVfdrQscSBUyZcdLUscCFVyIFiWOBCq1JATOfnTjpYNDoQqSSrdRs0s7RwICRiW941WLRscCAkY1pD3KYNlggMhAU0NOY9DsExwICSgyXMzWEY4EBLgyVosKxwICSicMjgQLP0cCAlwC8GywoGQAI9DsKxwICTA4xAsKxwICWhqyPuyo2WCAyEB7kOwrHAgJMCBYFnhQEhAYz7Htu6+89+apY8DIQH5HPSEA8HSz4GQgLxEd48DwdKvokCQdLKkZZKWS7pyB/vMkvSMpKWS/l+yZda3fC7nQLBMGHCyV0l54BbgRKADeFrSgoh4sWyfvYDvAydHxJuSxg1VwfUon8OBYJlQSQvhSGB5RLweEduAu4Ez+uxzDnB/RLwJEBHvJFtmfcvlRLf7ECwDKgmE8cCqsuWO4rpyBwF7S1osaYmk2f09kaQ5ktoltXd2dg6u4jqUl+hxC8EyoJJAUD/r+v71NwCfBk4DvgD8N0kHfeRBEfMjoi0i2saOHbvLxdarvFsIlhED9iFQaBFMLFueAKzuZ581EbER2CjpV8AM4JVEqqxz+ZyIgJ6eIJfrLz/N0qGSFsLTwBRJkyUNA84CFvTZ56fA5yU1SBoBHAW8lGyp9SuvQgi4lWBpN2ALISK6JF0GPAzkgdsiYqmki4vb50XES5IeAp4DeoAfRMQLQ1l4PeltFXT3BI35GhdjVoVKThmIiEXAoj7r5vVZ/kfgH5MrLT3yxUDwaEVLO49UTEDplMFXGizlHAgJ6D1l6PEHHi3lHAgJaMi5U9GywYGQgN4WQpebCJZyDoQE9PYhOA8s7RwICcgXX0WfMljaORASkCu1EBwIlm4OhATkc77saNngQEhA3lcZLCMcCAlwC8GywoGQAI9UtKxwICQg5xaCZYQDIQGlcQjuQ7CUcyAkwH0IlhUOhAT448+WFQ6EBPQGQpenc7OUcyAkIOdbqFlGOBASkPf9ECwjHAgJ8IebLCscCAnwh5ssKxwICWjIFV5GX3a0tHMgJKCYB3Q5ECzlHAgJ8DgEywoHQgL84SbLCgdCAnJuIVhGOBAS4BaCZYUDIQH+cJNlRUWBIOlkScskLZd05U72O0JSt6QvJ1di/XMgWFYMGAiS8sAtwCnAVOBsSVN3sN8/UJglerfieypaVlTSQjgSWB4Rr0fENuBu4Ix+9vsvwH3AOwnWlwoeqWhZUUkgjAdWlS13FNeVSBoP/AXwoSni+5I0R1K7pPbOzs5drbVu+ZTBsqKSQFA/6/r+5f8zcEVEdO/siSJifkS0RUTb2LFjK62x7pWuMjgPLOUaKtinA5hYtjwBWN1nnzbgbhXeGC3AqZK6IuKBRKqsc/l8bwvBn3+2dKskEJ4GpkiaDLwFnAWcU75DREzu/VnS7cDC3SUMoHwcQo0LMavSgIEQEV2SLqNw9SAP3BYRSyVdXNy+036D3UHvh5s8UtHSrpIWAhGxCFjUZ12/QRARF1RfVrp4pKJlhUcqJsBXGSwrHAgJkITkUwZLPwdCQhpy8g1SLPUcCAnJSR6paKnnQEhIPif3IVjqORASkpf84SZLPQdCQnI5nzJY+jkQEpLPuYVg6edASIj7ECwLHAgJycuBYOnnQEhIoYVQ6yrMquNASEgu55GKln4OhIT4lMGywIGQkJyvMlgGOBAS0pAT3b6HmqWcAyEhOY9UtAxwICQk75GKlgEOhIR4pKJlgQMhITlfZbAMcCAkJJ+TxyFY6jkQEpLPiS5fZbCUcyAkJC+3ECz9HAgJ8acdLQscCAkpjFSsdRVm1XEgJCQvTwdv6edASIhPGSwLHAgJcSBYFlQUCJJOlrRM0nJJV/az/VxJzxW/npQ0I/lS65tHKloWDBgIkvLALcApwFTgbElT++y2AviziJgOfA+Yn3Sh9c4TtVgWVNJCOBJYHhGvR8Q24G7gjPIdIuLJiPiP4uJvgAnJlln/3EKwLKgkEMYDq8qWO4rrduQi4N/62yBpjqR2Se2dnZ2VV5kCvmOSZUElgaB+1vX7ly/pOAqBcEV/2yNifkS0RUTb2LFjK68yBTxRi2VBQwX7dAATy5YnAKv77iRpOvAD4JSIWJtMeenh2Z8tCyppITwNTJE0WdIw4CxgQfkOkiYB9wPnR8QryZdZ/3L+tKNlwIAthIjoknQZ8DCQB26LiKWSLi5unwdcDYwBvi8JoCsi2oau7PrjFoJlQSWnDETEImBRn3Xzyn7+K+Cvki0tXYblc2zr8kwtlm4eqZiQ5sY8Wx0IlnIOhIQ0NeTo7gm6PJ+bpZgDISFNjYWXcotbCZZiDoSENDXkAdi6vbvGlZgNngMhIc3FFoL7ESzNHAgJKbUQHAiWYg6EhDQ1FPsQfMpgKeZASEiTTxksAxwICWl2p6JlgAMhIW4hWBY4EBLS26noPgRLMwdCQno7Fd1CsDRzICSkudGXHS39HAgJ+aCF4FMGSy8HQkI+6ENwC8HSy4GQkA+uMriFYOnlQEhI6ZTBLQRLMQdCQiQxrCHnTkVLNQdCgpoach6HYKnmQEiQb6NmaedASFBTQ86dipZqDoQENbkPwVLOgZCgpoa8P+1oqeZASFBzo1sIlm4OhAQVWggOBEsvB0KCmhrdqWjp5kBIUGEcglsIll4OhAQVxiG4hWDpVVEgSDpZ0jJJyyVd2c92SbqpuP05SYcnX2r982VHS7sBZ3+WlAduAU4EOoCnJS2IiBfLdjsFmFL8Ogq4tfh9t9LUkGf9li7eem8zI4fla12O7Yb2HN6IpEE/vpLp4I8ElkfE6wCS7gbOAMoD4QzgXyMigN9I2kvSfhHxh0FXlkKfPbCFO3/7Jp+9/tFal2K7qRV/f2pVj68kEMYDq8qWO/jov/797TMe+FAgSJoDzAGYNGnSrtZa904+9E94+NvH8sSrnUStizEbhEoCob/2R9+/90r2ISLmA/MB2traMvmeOXDcHhw4bo9al2E2KJV0KnYAE8uWJwCrB7GPmdW5SgLhaWCKpMmShgFnAQv67LMAmF282nA0sG536z8wy4IBTxkiokvSZcDDQB64LSKWSrq4uH0esAg4FVgObAIuHLqSzWyoVNKHQEQsovCmL183r+znAL6RbGlm9nHzSEUzK3EgmFmJA8HMShwIZlbiQDCzEgeCmZU4EMysxIFgZiUOBDMrcSCYWYkDwcxKHAhmVqLC55JqcGCpE3ijgl1bgDVDXE616r3Geq8PXGNSKqnxExExtr8NNQuESklqj4i2WtexM/VeY73XB64xKdXW6FMGMytxIJhZSRoCYX6tC6hAvddY7/WBa0xKVTXWfR+CmX180tBCMLOPiQPBzErqNhAGmmC2FiRNlPSYpJckLZX0reL6fST9QtKrxe9717jOvKTfS1pYj/UVa9pL0r2SXi6+nsfUU52S/rr4//gFSXdJaq51fZJuk/SOpBfK1u2wJklXFd8/yyR9oZJj1GUglE0wewowFThb0tTaVgVAF/A3EXEIcDTwjWJdVwK/jIgpwC+Ly7X0LeClsuV6qw/gX4CHIuKTwAwK9dZFnZLGA98E2iLiUArTD5xVB/XdDpzcZ12/NRX/Ls8CPlV8zPeL76udi4i6+wKOAR4uW74KuKrWdfVT508pzIq9DNivuG4/YFkNa5pQ/MM4HlhYXFc39RVrGA2soNipXba+Lurkg7lK96EwVcFC4KR6qA9oBV4Y6DXr+56hMK/KMQM9f122ENjx5LF1Q1IrcBjw78C+UZypqvh9XO0q45+BvwV6ytbVU30ABwCdwA+LpzY/kDSSOqkzIt4CbgDepDBh8bqI+Hm91NfHjmoa1HuoXgOhoslja0XSHsB9wLcj4v1a19NL0unAOxGxpNa1DKABOBy4NSIOAzZSH6cxABTPw88AJgP7AyMlnVfbqnbZoN5D9RoIdTt5rKRGCmHw44i4v7j6bUn7FbfvB7xTo/I+C3xR0krgbuB4Sf+njurr1QF0RMS/F5fvpRAQ9VLnfwJWRERnRGwH7gc+U0f1ldtRTYN6D9VrIFQywezHTpKA/w28FBH/VLZpAfDV4s9fpdC38LGLiKsiYkJEtFJ4zR6NiPPqpb5eEfFHYJWkg4urTgBepH7qfBM4WtKI4v/zEyh0etZLfeV2VNMC4CxJTZImA1OA3w74bLXotKmw8+RU4BXgNeDval1PsabPUWh2PQc8U/w6FRhDoSPv1eL3feqg1ll80KlYj/XNBNqLr+UDwN71VCfwXeBl4AXgDqCp1vUBd1Ho09hOoQVw0c5qAv6u+P5ZBpxSyTE8dNnMSur1lMHMasCBYGYlDgQzK3EgmFmJA8HMShwIZlbiQDCzkv8PCcNkFBFhFbIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from fipy import *\n",
    "\n",
    "u = 1.e-3\n",
    "L = 100.\n",
    "nx = 200\n",
    "dt = 200.\n",
    "dx = L/nx\n",
    "muo = 0.002\n",
    "muw = 0.001\n",
    "mesh = Grid1D(dx = L/nx, nx = nx)\n",
    "x = mesh.cellCenters\n",
    "sw = CellVariable(mesh=mesh, name=\"saturation\", hasOld=True, value = 0.)\n",
    "sw.setValue(1,where = x<=dx)\n",
    "sw.constrain(1.,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n",
    "\n",
    "eq = TransientTerm(coeff=1) + UpwindConvectionTerm(coeff = u\n",
    "*(sw**2./muw)/(sw**2./muw+(1-sw)**2./muo) * [[1]]) == 0\n",
    "\n",
    "sw.constrain(1.,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n",
    "\n",
    "steps = 100\n",
    "viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)\n",
    "for step in range(steps):\n",
    "    sw.updateOld()\n",
    "    swres = 1.0e6\n",
    "    while swres > 1e-5:\n",
    "        swres = eq.sweep(dt = dt, var = sw)\n",
    "        print(swres)\n",
    "        \n",
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Updates solution on October 8th, 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fipy import *\n",
    "\n",
    "# relperm parameters\n",
    "swc = 0.1\n",
    "sor = 0.1\n",
    "krw0 = 0.3\n",
    "kro0 = 1.0\n",
    "nw = 2.0\n",
    "no = 2.0\n",
    "\n",
    "# domain and boundaries\n",
    "u = 1.e-3\n",
    "L = 100.\n",
    "nx = 50\n",
    "dt = 200.\n",
    "dx = L/nx\n",
    "\n",
    "# fluid properties\n",
    "muo = 0.002\n",
    "muw = 0.001\n",
    "\n",
    "# define the fractional flow functions\n",
    "def krw(sw):\n",
    "    res = krw0*((sw-swc)/(1-swc-sor))**nw\n",
    "    return res\n",
    "\n",
    "def kro(sw):\n",
    "    res = kro0*((1-sw-sor)/(1-swc-sor))**no\n",
    "    return res\n",
    "\n",
    "def fw(sw):\n",
    "    res = krw(sw)/muw/(krw(sw)/muw+kro(sw)/muo)\n",
    "    return res\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "sw_plot = np.linspace(swc, 1-sor, 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the relative permeability and fractional flow curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3yV5f3/8dcne+9Fwggj7CEQUBFURAVxIIoVd9VKtbXDb4d+tbZ+21+tbb/229o6iqNqHRT3AATrRClK2ENG2CEJZJO9zvX74zrBCIGcQJL7nJPP8/E4j7Nuct6G+ObOdd/3dYkxBqWUUr4vwOkASimlOocWulJK+QktdKWU8hNa6Eop5Se00JVSyk8EOfXBSUlJJjMz06mPV0opn7R69epiY0xyW+85VuiZmZnk5OQ49fFKKeWTRGTv8d7TIRellPITWuhKKeUntNCVUspPaKErpZSf0EJXSik/0W6hi8gzInJIRDYd530RkUdEJFdENojIuM6PqZRSqj2e7KE/C8w4wfsXAVnu2zzg8VOPpZRSqqPaLXRjzKdA6Qk2mQU8b6yVQJyI9OqsgMeoPAhL7oGmhi77CKWU8kWdMYaeAexv9TzP/doxRGSeiOSISE5RUdHJfdr+lfDF4/D+/Sf355VSyk91RqFLG6+1uWqGMWa+MSbbGJOdnNzmlavtGz4LzvgefPEEbHz15L6GUkr5oc4o9DygT6vnvYH8Tvi6x3fBr6HvmfD2D+Dgli79KKWU8hWdUehvAze6z3Y5A6gwxhR0wtc9vsBguOpZCImChTdA3eEu/TillPIFnpy2+DLwH2CIiOSJyK0icruI3O7eZDGwC8gFngS+12VpW4tOs6Veuhve+h7o2qhKqR6u3dkWjTHXtPO+Ab7faYk6IvMsO/yy7D5Y8Qic9SNHYiillDfw/StFz/y+PVD67wdg93Kn0yillGN8v9BFYNajkDgIXr0ZDnft8VillPJWvl/oAKHRcPUL0FADC2+EpnqnEymlVLfzj0IHSB4Csx+HvFWw5G6n0yilVLfzn0IHO5Y++S5Y/Q9Y/ZzTaZRSqlv5V6EDnHc/DDwPFv8U8nTNUqVUz+F/hR4QCFc+DdG94F83QNUhpxMppVS38L9CB4hIgLkvQm0ZLLwJmhudTqSUUl3OPwsdIG0UXPZX2LcClv3C6TRKKdXl2r1S1KeNvgry18LKRyF9LIyZ63QipZTqMv67h97igl9D5hR4+4dwYI3TaZRSqsv4f6EHBsFVz0FUKiy4zq54pJRSfsj/Cx0gMhGueQnqyu10u3olqVLKD/WMQgd7kPTyx2D/F/YcdZ1uVynlZ/z7oOjRRsyGwk2w/H8hbTRMvM3pREop1Wl6zh56i6n3weAZ8N49Ot2uUsqv9LxCDwiAK+ZDwgB45SYo2+t0IqWU6hQ9r9ABwmJh7svQ3AQLroX6KqcTKaXUKeuZhQ6QNAjmPAOHtsAb3wWXy+lESil1SnpuoQNknQ8X/ha2vgsfP+h0GqWUOiU96yyXtpxxh91L//SPkDwURs1xOpFSSp2Unr2HDnZN0ov/BH0nwVvfh7zVTidSSqmTooUOEBQCV/8TolLsQVJdaFop5YO00FtEJsE1C6ChypZ6Q43TiZRSqkO00FtLHQFXPAn56+zwi04PoJTyIVroRxs6E87/FWx+HT7+ndNplFLKY3qWS1vO+jEU58Inv4fEQTD6W04nUkqpdukeeltE4JL/g36T7dDLvpVOJ1JKqXZpoR9Py5kvsX3sQdLS3U4nUkqpE9JCP5GIBLh2Ibia4aWrobbc6URKKXVcWujtSRpk99RLd8Ir34bmRqcTKaVUmzwqdBGZISLbRCRXRO5p4/1YEXlHRNaLyGYRubnzozqo/9lwyZ9h10ew5Od6OqNSyiu1e5aLiAQCjwIXAHnAKhF52xizpdVm3we2GGMuFZFkYJuIvGiMaeiS1E4YdwOU7IDP/2LnUp/0A6cTKaXUN3iyhz4RyDXG7HIX9AJg1lHbGCBaRASIAkqBpk5N6g2mPQDDZ8Gy+2HLW06nUUqpb/Ck0DOA/a2e57lfa+1vwDAgH9gI/MgYc8wE4yIyT0RyRCSnqKjoJCM7KCAAZv8demfD6/Ng/yqnEyml1BGeFLq08drRg8jTgXVAOnAa8DcRiTnmDxkz3xiTbYzJTk5O7nBYrxAcbud8iU6Dl+fq6YxKKa/hSaHnAX1aPe+N3RNv7WbgdWPlAruBoZ0T0QtFJsF1r4FphhevgppSpxMppZRHhb4KyBKR/iISAswF3j5qm33ANAARSQWGALs6M6jXSRoEc1+C8r3wr+uhqd7pREqpHq7dQjfGNAF3AkuBr4CFxpjNInK7iNzu3uw3wCQR2Qh8ANxtjCnuqtBeo98kuPxx2Pu5zs6olHKcR5NzGWMWA4uPeu2JVo/zgQs7N5qPGDUHyvfBB/8DcX1h2i+dTqSU6qF0tsXOMPkuO/Sy/GGISYcJ33E6kVKqB9JC7wwiMPNhqDwIi38G0b1g6MVOp1JK9TA6l0tnCQyCOU9D+lh49RbY/6XTiZRSPYwWemcKibSzM8ak29kZi3OdTqSU6kG00DtbZBJc/xpIALxwhR2GUUqpbqCF3hUSBsB1C6G6CF66CuornU6klOoBtNC7SsZ4uOo5KNwEC2+EJv+ZeFIp5Z200LvS4Avhskdg54fw5h3gOma+MqWU6jR62mJXG3s9VB2yFx5FJsOM39nTHJVSqpNpoXeHyXfZUv/icYhKgSn/5XQipZQf0kLvDiIw/UGoKf56T33cDU6nUkr5GS307hIQALMes1PtvvNDiEiEoTOdTqWU8iN6ULQ7BYXAt56HXqfBqzfD3hVOJ1JK+REt9O4WGgXXvQKxveGluVC40elESik/oYXuhMgkuOFNW+7/vAJKdjqdSCnlB7TQnRLXx5a6aYbnL4eKA04nUkr5OC10JyUPtvO+1JbBP2dDdYnTiZRSPkwL3WnpY+HaBXaBjBev1HlflFInTQvdG2ROtvO+FGyAl6+BxjqnEymlfJAWurcYMgNmPwF7PrMLZDQ3Op1IKeVjtNC9yehvwcw/wrZF7sm8mp1OpJTyIXqlqLeZeBs0VMG/H4DgcLj0EZ3MSynlES10bzT5LmiogU//AMERMOMhLXWlVLu00L3V1HuhoRpWPmrXKp32S6cTKaW8nBa6txKB6b+FxhpY/rDdUz/7p06nUkp5MS10byYCF//JlvqHv7F76mfc4XQqpZSX0kL3di3T7jbWwHv3QFAoZN/idCqllBfS0xZ9QWAQXPkMZE2Hd++CNf90OpFSygtpofuKlrnUB54Hb/8A1r3sdCKllJfRQvclwWEw9yXofza89T3Y8IrTiZRSXsSjQheRGSKyTURyReSe42xzroisE5HNIvJJ58ZURwSHwzULoO8keOO7sPkNpxMppbxEu4UuIoHAo8BFwHDgGhEZftQ2ccBjwGXGmBHAVV2QVbUIiYBr/wW9J8Br34Gv3nU6kVLKC3iyhz4RyDXG7DLGNAALgFlHbXMt8LoxZh+AMeZQ58ZUx2hZyi59LLzybdi2xOlESimHeVLoGcD+Vs/z3K+1NhiIF5GPRWS1iNzY1hcSkXkikiMiOUVFRSeXWH0tLMYukJE2Cv51A2xd7HQipZSDPCn0tiYRMUc9DwLGAxcD04H7RWTwMX/ImPnGmGxjTHZycnKHw6o2hMXCDW/YUl94I2xd5HQipZRDPCn0PKBPq+e9gfw2tnnPGFNtjCkGPgXGdE5E1a7wOLjxTeg12pb6V+84nUgp5QBPCn0VkCUi/UUkBJgLvH3UNm8BU0QkSEQigNOBrzo3qjqhlj31XqfZMfUtR/8VKaX8XbuFboxpAu4ElmJLeqExZrOI3C4it7u3+Qp4D9gAfAk8ZYzZ1HWxVZtaSj19HLx6M2x5y+lESqluJMYcPRzePbKzs01OTo4jn+336g7Di3MgLwfmPA0jZjudSCnVSURktTEmu6339EpRf9Ry9kvvCfDqrXpFqVI9hBa6vwqNtqXebxK8fptO6KVUD6CF7s9Co+DahTBwKrx9J6x6yulESqkupIXu70IiYO7LMPgiWPQT+M+jTidSSnURLfSeIDjMTr077DJYeq9d0k4p5Xe00HuKoBCY8w8YdRV88Gv48Lfg0BlOSqmuoUvQ9SSBQTD773YZu0//AA1VMP1Bu3apUsrnaaH3NAGBcOlfISQKVj4G9Yfh0kfs60opn6aF3hMFBMCMhyA0xu6p11fBFU/aYRmllM/SQu+pROC8++xFSMt+YYdfvvVPe1aMUson6UHRnm7SD+yQS+4H8MKVUFfhdCKl1EnSQlcw/iY750vel/DcpVCli48o5Yu00JU18kqY+xIUbYN/zICyvU4nUkp1kBa6+trg6XDDm1BdBM9Mh4NbnE6klOoALXT1Tf3OhJuX2IuO/jED9n3hdCKllIe00NWxUkfArcsgIgmenwXblzmdSCnlAS101bb4fnDLUkgeDC/PhfULnE6klGqHFro6vqhkuOldyDwL3vgufP6Izv+ilBfTQlcnFhYD171ql7F7/3547x5wNTudSinVBr1SVLUvKBSufAai02Hlo1BZALPn22l5lVJeQwtdeSYgAGY8CDHpsOw+e/HR3BchIsHpZEopNx1yUR0z6U6Y8wwcyIFnZkD5fqcTKaXctNBVx428Eq5/HSoL4anzoXCj04mUUmihq5PVfwrc8p6dR/2ZGbDjfacTKdXjaaGrk5c6HL7zb0gYAC9dDauedjqRUj2aFro6NTHpdqqAQefDov+CpfeBy+V0KqV6JC10depCo+xMjRNug//8DRbeAA01TqdSqsfRQledIzAIZv4Rpv8Oti6C5y6BqkNOp1KqR9FCV51HBM78Hlz9gp1698nzoHCT06mU6jG00FXnG3YJ3LwYXE12XvVtS5xOpFSPoIWuukbGOLjtQ0gcBC9fA5//RSf2UqqLeVToIjJDRLaJSK6I3HOC7SaISLOIzOm8iMpntZwBM3wWvP9LeOv70FTvdCql/Fa7hS4igcCjwEXAcOAaERl+nO1+Dyzt7JDKh4VEwJx/wDn3wLoX7YIZ1cVOp1LKL3myhz4RyDXG7DLGNAALgFltbPcD4DVAT21Q3xQQAFP/G658GvLXwpNTdboApbqAJ4WeAbSegSnP/doRIpIBzAaeONEXEpF5IpIjIjlFRUUdzap83ag59mBpcyM8fSFsfsPpREr5FU8KXdp47eijW38G7jbGnHDlA2PMfGNMtjEmOzk52dOMyp9kjId5H0PqSHjl2/DBr3XBDKU6iSfzoecBfVo97w3kH7VNNrBARACSgJki0mSMebNTUir/Ep0G334XFv8Mlj9sz1W/8kkIi3U6mVI+zZM99FVAloj0F5EQYC7wdusNjDH9jTGZxphM4FXge1rm6oSCQuHSv8DFD8POD+DJaVC8w+lUSvm0dgvdGNME3Ik9e+UrYKExZrOI3C4it3d1QOXHRGDCd+DGt6G2zF5ZunWx06mU8lliHLrYIzs72+Tk5Djy2coLle+3k3rlr4UpP4Wp99q51pVS3yAiq40x2W29p1eKKu8Q1wdufg/G3QjL/xdeuBKqS5xOpZRP0UJX3iM4DC77K1z6COxdAfPPgQNrnE6llM/QQlfeZ/xNdnk7sJN7rX7O2TxK+QgtdOWdMsbBvE+g31nwzg/hze/pohlKtUMLXXmvyES4/jU4++ew7iV4ahoUbXc6lVJeSwtdebeAQDjvPlvsVYdg/rmwYaHTqZTySlroyjcMmga3L4deY+D12+CdH0FjrdOplPIqWujKd8Skw03vwOS7YPWz8NQFULLT6VRKeQ0tdOVbAoPg/Afg2lfgcB78/WxYv8DpVEp5BS105ZsGXwi3f2aHYN74Lrw+D+ornU6llKO00JXviu1th2DO/W/Y+IrdW89f63QqpRyjha58W0AgnHsPfHuRXa/0qQtgxV/B5XI6mVLdTgtd+Yd+k+wQzODpsOwX8OIcqCx0OpVSxzh4uI5DlXVd8rW10JX/iEiAq1+wc6zv/RweOxO+esfpVErR2Oxi6eZCbn12FZMe+pCnl+/uks/xZMUipXxHyxzrmVPs+er/uh7G3gAzHoLQKKfTqR4m91AVC3P28/qaPIqrGkiJDuW7Zw/gW9l92v/DJ0ELXfmn5CFw67/h49/BZ/8Hez6DK+ZDn4lOJ1N+rrKukUUbCnhldR6r95YRFCCcNzSFqyf04ZzByQQFdt3AiBa68l9BIXD+ryDrAnj9u3bmxrN/Zm+BwU6nU37E5TKs3FXCK6vzWLKpgLpGFwOTI7l35lBmj+1NcnRot+TQQlf+r98kuOMzWHI3fPJ72P4eXP4EpA53OpnycftKanh1TR6vrc7jQHkt0WFBXDGuN1eN781pfeIQkW7No4WueoawWJj9BAy9GN75sV0849x7YNKP7NWnSnnocF0jizcU8NqaPFbtKUMEJg9K4uczhjB9RBphwc4tnag/yapnGXYp9D0TFv0XfPBruyj15Y9D8mCnkykv1uwyLN9RxGtrDrBscyH1TS4GJEfys+lDmD02g/S4cKcjAlroqieKTIKrnoNNr8Hin8Lfp8B598MZd+jC1OoIYwyb8w/z5toDvLU+n6LKemLDg/lWdh+uGJfhyJBKe7TQVc8kAqPmQOZkOwSz7D7Y8hbM+ps9Q0b1WPnltby1Lp831uax/WAVwYHCuUNSuGJsBucNSyE0yHv/0RdjjCMfnJ2dbXJychz5bKW+wRi7aMZ7d0NDNZzzczjrx3omTA9SUdPIkk0FvLUun5W7SzAGxvWNY/a43lwyqhfxkSFORzxCRFYbY7Lbek/30JUSgTFXw8DzYMnP4MP/B5vftHvr6WOdTqe6SF1jMx9uPcSbaw/w8bYiGppd9E+K5EfTspg9NoN+iZFOR+wwLXSlWkQlw1XPwsg5sOgn8OQ0mHSnnc0x2DsOeqlT09TsYsXOEt5en8/STYVU1jeRHB3K9Wf04/Kx6YzKiPW6cfGO0EJX6mjDLrFj6+/fD5//xY6tX/wnuwye8jkulyFnbxnvrM9n8cYCSqobiAoNYsbINC4/LYMzByYSGOC7Jd6aFrpSbQmPg8v+CqOugnfvgheusHvu0x+E6FSn06l2GGPYeKCCdzcU8O76fPIr6ggLDmDasFQuHZ3OuUOSHT1fvKtooSt1Iv3PhjtW2Plglj8MO9630wmMvxkCdLJSb9JymuGijQUs2lDAvtIaggOFs7OSufuioUwblkpUqH9Xnp7lopSninfYvfU9y6H3BLjkz5A20ulUPZoxhq2FlSzaUMCijQXsLq4mMEA4a1ASl4zuxfThacRG+NfZSnqWi1KdISnLLnm34V+w9F675N3E2+xB0/A4p9P1GC174os3FrBkUyG7i6sJEJg0MIl5Zw9g+og0ErzoNMPu5FGhi8gM4C9AIPCUMeaho96/Drjb/bQKuMMYs74zgyrlFURgzFzIutCe3vjF3+0Vp+f/D4y5RodhukjLmPjijYUs2VTA3pIaAgOEMwckctuUAVw4IpWkqO6Z0dCbtTvkIiKBwHbgAiAPWAVcY4zZ0mqbScBXxpgyEbkIeMAYc/qJvq4OuSi/kL/OTh+Qtwp6T4SZf4T005xO5ReaXYbVe8tYsqmApZsKya+oIyhAmDQoiZkj07iwh+6Jn+qQy0Qg1xizy/3FFgCzgCOFboxZ0Wr7lUDvk4+rlA9JPw1uWQbrX4L3fwXzz4Xsm2HqLyAy0el0PqehycXKXSW8t7mQZZsPUlxVT0hQAGdnJXHXBYO5YHgqcRE9r8Q95UmhZwD7Wz3PA060930rsKStN0RkHjAPoG/fvh5GVMrLBQTA2Oth6CXw0YOw6inY+BqcezdMuM0utKGOq7q+iU+2F7F0cyEfbj1EZV0TESGBTB2awowRaUwdmuL3Z6d0Fk++S22dcd/mOI2ITMUW+uS23jfGzAfmgx1y8TCjUr4hPA5m/gGyb7EHTZfeC6uehum/hcEz7Pi7AqCosp4Ptx5k6eaDfJZbTEOTi/iIYGaMSGP6iDQmZyX55XniXc2TQs8DWq9o2hvIP3ojERkNPAVcZIwp6Zx4SvmglKFww+v2nPWl98LLc6H/OfaipB56mqMxhp1FVSzbcpB/bznI2v3lGAMZceFcf3o/LhyRSna/+C5db7Mn8KTQVwFZItIfOADMBa5tvYGI9AVeB24wxmzv9JRK+aKsC2DAuZDzD/j4QTvv+phrYeq9EJvhdLou19jsYvXeMj746iDvbznInpIaAEZlxHLX+YOZNiyF4b1ifHruFG/TbqEbY5pE5E5gKfa0xWeMMZtF5Hb3+08AvwQSgcfcfzlNxzsKq1SPEhgMp8+zc68vfxi+nA+bXoXTb4fJd/nd+etl1Q18sr2ID7Ye4pNthzhc10RIYABnDkzkO1MGMG1YCr1idaKzrqJXiirVncr22vPXNy6E8HiY8lN7cVKQb55DbYxh28FKPtpaxEdbD5GztxSXgaSoEKYOSWHasFQmZyXpQc1OdKLTFrXQlXJCwXp7muOujyC2r12wevTVPrFgdU1DE5/nlvDRtkN8vPUQ+RV1AAzrFcP5w2yJj86IJcBPZjD0NlroSnmrnR/Cvx+wBZ+YBVP/G4bP9qorTlsOaH68rYhPthfxxa5SGppdRIYEMjkrialDUjh3SAppsWFOR+0RtNCV8mbGwFfvwEe/haKtkDoSpt4HQy5y7FTHyrpGVuws4ZPtRXyyrYgD5bUADEqJ4pzByZw3NIUJmQmEBHnPPzw9hRa6Ur7A1WznhfnoQSjbDRnj7RkxA6d1ebE3u+xcKcu3F7F8RzFr9pXR5DJEhgRy1qAkzh2SwtmDk+gdH9GlOVT7tNCV8iXNjbD+ZfjkD1CxH9LHwTl3w+DpnVrsB8pr+WxHEZ/uKObz3GLKaxoBGJkRw9lZyUzJSmZ8v3jdC/cyOn2uUr4kMBjG3Qij59o5Ypb/CV6+GtJGwzk/hyEXn9QYe0VNI//ZVcxnucV8nlvC7uJqAFJjQjl/WCpTspKYPCiJRJ210GfpHrpS3q65ETYstOexl+6ElOEw5Scw/PITnhVT29BMzt5SVuwsYUVuMRsPVOAyEBESyBkDEjlrkC3wwalRenGPD9EhF6X8QXMTbH4DPv0jFG+D+Ew480447ToIiaChycX6vHJW5JawYmcxa/eV09DsIihAGNMnjsmDkpiclcSY3nE6jOLDtNCV8icuF2xbjOuzPxNwYBW1wXEsCr+Mh8vOpqAxAhEYkR7DpIFJnDkwkQmZCXphjx/RMXSl/EB9UzMb8ir4cncpK3clk7PvJ4xq2sx3m99lTuPzXBa8kPwhV5Ew9QfEZAxxOq5ygBa6Ul6qpqGJNXvL+XJ3CV/sLmXt/nIamlwADEmN5uoJfTljwGmM7f9DqNpByOePkLnpZch9AYbMhDPugMzJOm1vD6JDLkp5iaLKelbvLWXVnjJy9pSyKf8wzS5DgMDIjFgmZiYwsX8CEzITiD/e0muHC+wCGznPQG0ppI6yxT7ySgjWKzn9gY6hK+VlXC7DruIqVu8tI2dPGTl7y46cRhgaFMCYPnFk94tnYv8ExveLJzosuGMf0Fhrz4xZ+TgUfQWRyTD2Brs8XpyuFubLtNCVclh1fRPr88pZs7eM1XvLWLOvnIpaeyFPfEQw4/slMCEznuzMBEZmxBAa1Emr9RgDuz6GL/4OO5ba54OnQ/atMGgaBOiqQL5GD4oq1Y1a9r7X7Ctn7b5y1u4rY/vBSlzufafBqVHMHJXGuL7xjO8XT/+kyK47D1wEBk61t/L9sPpZWPM8bH/P7qmPv9muhxqV0jWfr7qV7qErdYoOVdaxfn8F6/eXsz6vnHX7y6msawIgOiyI0/rEMbZvPGP7xjGuTzyxER0cPulsTQ2w9V07zr5nOQQEQdZ0GHcDDLrAJ6bw7cl0D12pTlJR08im/Ao25FWwIa+c9fvLj8wHHhggDEmN5pLR6Yzra0t8QFKk980LHhQCI6+wt6LtsPafdu6YbYsgKhXGXGP32pOynE6qOkj30JU6jsN1jWw+cJhNByrYcKCCjXnlR9bFBOibEMFpfeIY0yeOMb1jGZEeS3iIj45JNzfCjmWw9gXYvhRMM/SeCKO/BSNmQ2SS0wmVmx4UVaodpdUNbDpQwab8Clvi+RXsbVXeGXHhjMqIZVTvWEb3jmVURixxEcc5ddDXVRbC+gX2LJlDm+2QzMBpttyHzIQQnULXSVroSrm5XIZ9pTVsKTjMlvzDR+4LD9cd2aZPQjgj02MZmRHL8PQYRmXEktRTZyAs3GTXP934Khw+ACFRMPRiOzHYwPP03HYHaKGrHqmitpFthZVsLTzM1sJKthYcZlthJdUNzYAd8x6UHMXw9BiG9YpmZEYsI3rFOn/Q0hu5XLD3c9jwL7u6Ul05hETbVZVGXG734LXcu4UWuvJrdY3N5B6qYvvBSrYftPfbCiuPLJsGEBsezJC0aIalRTM8PYbhvWLJSo0iLNhHx7yd1NwIuz6BLW/as2Vqy+ye++AZMHQmDDofwmKdTum3tNCVX6htaGZnURU7i6rYcbCKHYcq2XGwij0l1UfO8Q4OFAYmRzEkLZqhaTEMTYtmaK9o0mLCdM7vrtDcCLs/dZf7IqgpsWPumZPtePvgGRDfz+mUfkULXfkMYwyl1Q3sLKpml7u8dxZVk3uoiv1lNbT8uAYGCP0SIxicEs3gtGiGpEYzJC2KfomRBAfqXN+OcDVDXg5sWwzbltg52wFSRkDW+XZYpu8ZENRDj0d0Ei105XVqGprYU1zDnpJqdhdXs6uomt3FVewqrj6ytiXYeU36J0UyKCWKrJRoe58aRWZipC7S4O1Kdtpi3/4e7FsJrkYIjrR774Om2YJPHKizQXaQFrpyRFV9E3tLqtlbUuO+VbPH/bygou4b26bGhJKZGMnAlCgGJkcxMDmSgclRpMeFE+htF+aojquvhN3LYecHkPsBlO22r8f0tgXffwpkTtHhGQ/olaKqSzQ2uyisqGN/aQ37SmvYX1bDvtJa9pfWsL+0hpLqhm9snxQVQr/ESM4ckEj/pEj6J0fSPymSzMRIInVFHf8WGh0LiMYAAAnVSURBVG0PmA6daZ+X7ISdH9qpB3Lfhw0L7OuxfW259z0T+pwOiYNOakHsnkr30NVxVdc3UVBRS355HfnltRworyWvrJYDZbXkldVQeLjuyMFIsOPaGXHh9EkIp29CBH0SIshMjKRfYgT9EiN1GTTVNmOgaKvdg9+zHPZ8ZudyBwiPh94ToM9Ee+VqxngIjXI2r8N0D10do6q+icKKWgor6imoqOXg4ToKKuoorKgjv8IWeMv0ri0CBHrFhpMRH84ZAxLpHW8f94m35d0rNowgPSCpOkoEUobZ2+nz7DnvJbmw/wvI+xL2f2mnJbAb2zlmep0G6afZ+16j7W8ASvfQ/YnLZSiraaC4qoHiqnqKKu3t4OE6DrnvW563XFzTWnxEMKkxYWTEhdMrLoz0uHDSY8PpFWsfp8WG6Rkkyhm1ZfYMmgOrIX8dFKyDygL3mwIJAyB1OKS0uiUM8MuZI3UP3Uc1uwwVtY2UVjdQXtNAabW9lVQ3UFLVQGl1/ZHHxVX2cbPr2H+gw4IDSI0JIyU6lGHpMZwzJJnUmDB6xYZ9414vslFeKzwesi6wtxaVB6FgvS33gvVw6Ct7Lryx664SGAJJg+04fOJASBj49eOIRL88u8ajQheRGcBfgEDgKWPMQ0e9L+73ZwI1wLeNMWs6OatPqm9qprKuyX1r5HCtva+sa+JwXSPlNY1U1H59K69t5LC7xA/XNXK8X6AiQgJJjAohITKUtNgwRmbEkBwdSlJU6Dfuk6NDiQ4N0otqlP+JToXoC2HwhV+/1lgLxdttuR/aYu8LN9rpCkyr30pDYyG+L8T2cd96Q5z7cUy6XbIv0PemgGi30EUkEHgUuADIA1aJyNvGmC2tNrsIyHLfTgced997JZfL0Ohy0dRsaGx20dDkot59a2hy0eB+ra6xmdrGZuoam6lvdFHX1ExtQzM1Dfb1moYmaurt85rGZqrrm6iut+Vd3WAfNzafeEgrQCAmPJjY8GDiwoOJCQ+mT3w4CZEhxEWEEB8R/I3HiVGhJEaG6N60Um0JDodeY+ytteZGKN9nz64pyYXSnXYFp7I99mBsQ+WxXys8wc4PH5Vi7yOTICwOwuPc9/Hux7EQHOG+hUNQmGNn5niyhz4RyDXG7AIQkQXALKB1oc8Cnjd2QH6liMSJSC9jTMGxX+7UfLztEL95dwsGwIDBXl1o7+0whcuYb9y33BpdhqZmF22MSnRYWHAAESFBhAcHEhFib1FhQSRGRhAVGkRUWBCRoUFEhQYRHea+hQa7H9v7mPBgokODvG8BBKX8TWCwHWpJHAhceOz7teVQkQcV++3YfNUh9+2gvd//hZ3WoKHKs88LCrflHhBk120NCAIJ+PrxuJtg0p2d+p8InhV6BrC/1fM8jt37bmubDOAbhS4i84B5AH37ntzK49FhwQxNiwEBsV/TfW+fBwQIgSIEBsg3HgcGCEGBQnBAgL0PDCAowN6HBAUQGnT0fSBhwfY+PCSQsOBAwoICCAsOJDw4UEtYKX8S7t7zTht54u2aG6Guwv4DUFf+9X1THTTUQGONHfZpuXc12aEeV8vN/byL1nD1pNDbaq6j93E92QZjzHxgPtizXDz47GOM72cX1lVKqW4XGGyHXrx0BSdPBnrygD6tnvcG8k9iG6WUUl3Ik0JfBWSJSH8RCQHmAm8ftc3bwI1inQFUdMX4uVJKqeNrd8jFGNMkIncCS7GnLT5jjNksIre7338CWIw9ZTEXe9rizV0XWSmlVFs8Og/dGLMYW9qtX3ui1WMDfL9zoymllOoIvY5bKaX8hBa6Ukr5CS10pZTyE1roSinlJxybPldEioC9J/nHk4DiTozTWbw1F3hvNs3VMZqrY/wxVz9jTHJbbzhW6KdCRHKONx+wk7w1F3hvNs3VMZqrY3paLh1yUUopP6GFrpRSfsJXC32+0wGOw1tzgfdm01wdo7k6pkfl8skxdKWUUsfy1T10pZRSR9FCV0opP+HVhS4iM0Rkm4jkisg9bbw/VET+IyL1IvJTL8p1nYhscN9WiMiYtr6OA7lmuTOtE5EcEZnsDblabTdBRJpFZI435BKRc0Wkwv39Wiciv/SGXK2yrRORzSLyiTfkEpGftfpebXL/XSZ4Qa5YEXlHRNa7v1/dMhusB7niReQN9/+TX4pIO8slecAY45U37FS9O4EBQAiwHhh+1DYpwATgt8BPvSjXJCDe/fgi4AsvyRXF18dNRgNbvSFXq+0+xM7qOccbcgHnAu92x89VB3PFYdf07et+nuINuY7a/lLgQ2/IBdwL/N79OBkoBUK8INcfgV+5Hw8FPjjVz/XmPfQji1MbYxqAlsWpjzDGHDLGrAIavSzXCmNMmfvpSuwKTt6Qq8q4f3qASNpYJtCJXG4/AF4DDnVDpo7k6m6e5LoWeN0Ysw/s/wdekqu1a4CXvSSXAaJFRLA7NaVAkxfkGg58AGCM2QpkikjqqXyoNxf68RaedlpHc90KLOnSRJZHuURktohsBRYBt3hDLhHJAGYDT9B9PP17PNP9q/oSERnhJbkGA/Ei8rGIrBaRG70kFwAiEgHMwP4D7Q25/gYMwy6LuRH4kTHG5QW51gNXAIjIRKAfp7jz582F7tHC0w7wOJeITMUW+t1dmsj9cW281tZC3W8YY4YClwO/6fJUnuX6M3C3Maa5G/K08CTXGuy8GWOAvwJvdnkqz3IFAeOBi4HpwP0iMtgLcrW4FPjcGFPahXlaeJJrOrAOSAdOA/4mIjFekOsh7D/M67C/oa7lFH9z8GjFIod468LTHuUSkdHAU8BFxpgSb8nVwhjzqYgMFJEkY0xXTl7kSa5sYIH9jZgkYKaINBljurJA281ljDnc6vFiEXnMS75feUCxMaYaqBaRT4ExwHaHc7WYS/cMt4BnuW4GHnIPN+aKyG7smPWXTuZy/3zdDOAeDtrtvp28rj5ocQoHFYKAXUB/vj6oMOI42z5A9x0UbTcX0Be7vuokb/p+AYP4+qDoOOBAy3Nv+Ht0b/8s3XNQ1JPvV1qr79dEYJ83fL+wwwcfuLeNADYBI53O5d4uFjtGHdnVf4cd+H49Djzgfpzq/rlP8oJccbgPzgK3Ac+f6ud67R668WBxahFJA3KAGMAlIj/GHkk+fNwv3A25gF8CicBj7r3OJtPFM755mOtK4EYRaQRqgauN+6fJ4VzdzsNcc4A7RKQJ+/2a6w3fL2PMVyLyHrABcAFPGWM2OZ3LvelsYJmxvz10OQ9z/QZ4VkQ2YodC7jZd+1uWp7mGAc+LSDP2rKVbT/Vz9dJ/pZTyE958UFQppVQHaKErpZSf0EJXSik/oYWulFJ+QgtdKaX8hBa6Ukr5CS10pZTyE/8fp83zLvH1rpgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXiU9b3+8fcnIQlbgMSEnQSQRXYIIaCeWq1aUY/FrRVQQVyoVjkuP1s9ntau52q17altXRAtKm6gFautVrRUxVpFEgirLCEIJCBJIBvZJpn5/v5ItCkGmcAkz8zkfl1XruSZeZK5HWZuv/N9NnPOISIikS/G6wAiIhIaKnQRkSihQhcRiRIqdBGRKKFCFxGJEp28euCUlBQ3ePBgrx5eRCQi5eTklDjnUlu6z7NCHzx4MNnZ2V49vIhIRDKz3Ue7T1MuIiJRQoUuIhIlVOgiIlFChS4iEiVU6CIiUeKYhW5mi82syMw2HeV+M7PfmVmemW0ws4zQxxQRkWMJZoT+JDD9S+4/Hxje9DUfeOTEY4mISGsdcz9059wqMxv8JavMAJa4xvPwfmhmvcysn3Nuf4gyioiEtUDAUVFbT2l1PaXVPsqr6ymr8eFrCODzO+obAtT7G798fkdmehJnjGjx2KATEooDiwYAe5stFzTd9oVCN7P5NI7iSUtLC8FDi4i0vboGP7sPVrP3UNNXaQ0FpdXsPVTD/vIaymrqac2lJW786slhW+jWwm0t/qc55xYBiwAyMzN1ZQ0RCTvlNfVs2VfBlv0VbN5XzpZ9FeQVHaYh8K/K6hwXw6CkrgxK7kpGei+Su8bTq2s8vbrGkdT0vWeXODrHxRIXG0N8bAydYo242BjiYg2zlmrzxIWi0AuAQc2WBwL7QvB3RUTaXMnhOj7YeZAP8g/ywc6D7Cqp+vy+3okJjO7fg7NH9WZEn0TSkrsyMKkrKd3j26yUT0QoCv1V4BYzWwpMBco1fy4i4aq23s97O0p4P6+ED3YeZNuBSgC6J3Ri6pBkLp88kDH9ezCmf09SExM8Tts6xyx0M3seOBNIMbMC4IdAHIBzbiHwOnABkAdUA/PaKqyIyPGo8fl5d3sRr238lL9/fIAqn5/OcTFMGZzMjEn9Oe3kFMb270Gn2Mg+NCeYvVxmHeN+B9wcskQiIiFQ1+Bn5cdFvLZxP29vLaLa5ye5WzzfmNif88f2Y+rQZBI6xXodM6Q8O32uiEhb2H2wiuc+2sOL2QUcqvKR0j2eSyYN4MJx/cgakhzxo/Avo0IXkYjX4A+wcmsRz3y4m/d2lBAbY5wzqjdXTk3n9GEpxMaE3wbMtqBCF5GIVVXXwDMf7uaJ9z/h04pa+vbozO3njOCKKYPo27Oz1/HanQpdRCJOZW09Sz7YzePv5VNaXc9pJ5/Ej2eM4exTekf1lMqxqNBFJGKU19Tz5Puf8Id/5FNR28BZI1NZcPZwMtKSvI4WFlToIhL2anx+Hnsvn8dW5VNZ18C5o/uw4GvDGD+wl9fRwooKXUTClnOO1zbu5+evb6WwrIbzxvTh1rNHMLp/D6+jhSUVuoiEpU2F5fzkz1v46JNDjO7Xg//71gSmDj3J61hhTYUuImGl5HAdv35zG0vX7CWpazw/v3Qc38oc1GF2PTwRKnQRCQvOOV5eV8iPXt1Mtc/PdacPYcHZw+nZJc7raBFDhS4iniuqqOWelzfyt4+LmJyexH2XjWdY7+5ex4o4KnQR8UzzUXldQ4DvXziKeacP0fTKcVKhi4gnDlTUcs/yjazc2jgq/+Xl4xmaqlH5iVChi0i7e2vLAe58cT219X6NykNIhS4i7abeH+BXK7bx6Kp8xg7owe9mTtKoPIRU6CLSLvaX13DLc+vI2V3KVdPS+P6Fo+kcF13nI/eaCl1E2ty724u5fVkutfV+fjtzIjMmDvA6UlRSoYtIm/EHHA/8bTsPvp3HiN6JPHRlhnZHbEMqdBFpE4frGrj1+XWs3FrENycP5CczxtIlXlMsbUmFLiIht6+shmufXMOOosP8ZMYY5pw62OtIHYIKXURCav3eMq5fkk2tz8/ia6bw1RGpXkfqMFToIhIyr2/cz+3LcklNTODZ66cyok+i15E6FBW6iJww5xwPv7OTX67YxuT0JB69ejIp3RO8jtXhqNBF5IT4A47v/2kjz3+0lxkT+3PfZeO1f7lHVOgictx8DQFuX5bLaxv3c/NZJ3Pn10dipkP4vaJCF5HjUu1r4MZn1rJqezH/c8EobjhjqNeROjwVuoi0WnlNPdc9uYa1e0q577JxXDElzetIggpdRFqpuLKOuYs/YkdRJQ/OzuCCcf28jiRNVOgiErTCshquenw1+8treHyu9jEPNyp0EQlKQWk1Mxd9SHlNPc9cN5XMwcleR5IjqNBF5Jj2l9cw+7HVlNfU89z10xg3sKfXkaQFMcGsZGbTzWybmeWZ2d0t3N/TzP5sZuvNbLOZzQt9VBHxwoGKWmYt+pDSKh9PXzdVZR7GjlnoZhYLPAScD4wGZpnZ6CNWuxnY4pybAJwJ/NrM4kOcVUTaWVFlLbMe+5DiyjqevDaLiYN6eR1JvkQwI/QsIM85l++c8wFLgRlHrOOARGs8oqA7cAhoCGlSEWlXJYfrmP3Yaj4tr+XJa7OYnJ7kdSQ5hmAKfQCwt9lyQdNtzT0IjAL2ARuBW51zgSP/kJnNN7NsM8suLi4+zsgi0tYOVfm46vHVFJRWs/iaKUzRBtCIEEyht3Qcrzti+TwgF+gPTAQeNLMeX/gl5xY55zKdc5mpqdrdSSQcVdTWc9Xjq9lVUsUf5k5h2tCTvI4kQQqm0AuAQc2WB9I4Em9uHrDcNcoDdgGnhCaiiLSX2no/1z+VzY6iSh69ejKnD0vxOpK0QjCFvgYYbmZDmjZ0zgRePWKdPcDZAGbWBxgJ5IcyqIi0rQZ/gP96fh0f7TrEr745gTNH9vY6krTSMfdDd841mNktwAogFljsnNtsZjc23b8Q+CnwpJltpHGK5i7nXEkb5haREHLO8T8vb+LNLQf44UWjmTHxyM1kEgmCOrDIOfc68PoRty1s9vM+4OuhjSYi7eWXK7axLHsvC742jHmnD/E6jhynoA4sEpHo9Yd/7OLhd3YyKyuNO84d4XUcOQEqdJEO7E/rCvnpX7YwfUxffnbxWF2cIsKp0EU6qH/sKOHOF9czbWgyD8ycSGyMyjzSqdBFOqDtByq56ZkcTk7tzqI5mboGaJRQoYt0MMWVdcx7Yg2d42NZPG8KPTrHeR1JQkSFLtKB1Pj8XL8km4NVdfxhbiYDenXxOpKEkM6HLtJBBAKO25flsqGgjEevmsz4gTpzYrTRCF2kg7jvja28sflTvn/haL4+pq/XcaQNqNBFOoBnV+/m0VX5zDk1nWtPH+x1HGkjKnSRKPfejmLufWUzZ41M5d7/HK19zaOYCl0kiu0qqeLmZ9cyvHd3fj87g06xestHM/3rikSpitp6bliSTWyM8dicTLonaB+IaKd/YZEo5A84bluayyclVTx93VQGJXf1OpK0A43QRaLQr97cxt+3FvHDb4zh1JN1xaGOQoUuEmVeyS3kkXd2MntqGldPS/c6jrQjFbpIFNlQUMb3/riBrCHJ/OiiMV7HkXamQheJEkUVtcxfkkNK9wQeuTKD+E56e3c02igqEgV8DQG+8+xaymvqeemm0zipe4LXkcQDKnSRKPC/r20he3cpD86exOj+PbyOIx7RZzKRCPdSTgFPfbCbG74yhP8c39/rOOIhFbpIBNtUWM49L29k2tBk7pp+itdxxGMqdJEIVVrl48ZnckjuFs+DOqxf0By6SETyBxy3LsulqKKOF248lRRtBBVU6CIR6TdvbWfV9mJ+fuk4Jg7ShSqkkT6jiUSYFZs/5cG385g5ZRCzstK8jiNhRIUuEkF2lVRx5wvrGT+wJz/6ho4ElX+nQheJEDU+Pzc9k0NsrPHwlRl0jov1OpKEGc2hi0SIe1/ZxLYDlTxxzRQGJul0uPJFGqGLRIBla/bwYk4BC84axpkje3sdR8KUCl0kzG0qLOcHr2zmK8NTuPWcEV7HkTAWVKGb2XQz22ZmeWZ291HWOdPMcs1ss5m9G9qYIh1TeU0933l2Lcld43ngionExugCz3J0x5xDN7NY4CHgXKAAWGNmrzrntjRbpxfwMDDdObfHzPSZUOQEOee488X17CurYdm3T9UZFOWYghmhZwF5zrl855wPWArMOGKd2cBy59weAOdcUWhjinQ8j67K560tB7jnglFMTk/yOo5EgGAKfQCwt9lyQdNtzY0AkszsHTPLMbM5Lf0hM5tvZtlmll1cXHx8iUU6gNX5B/nlim1cOK4f804f7HUciRDBFHpLk3buiOVOwGTgQuA84Adm9oWtN865Rc65TOdcZmpqaqvDinQEJYfrWPD8OtKSu/KLy8ZhpnlzCU4w+6EXAIOaLQ8E9rWwTolzrgqoMrNVwARge0hSinQQ/oDjtqW5lNfU8+S8LBI7x3kdSSJIMCP0NcBwMxtiZvHATODVI9Z5BfiKmXUys67AVODj0EYViX4P/j2Pf+SV8JMZY3TlIWm1Y47QnXMNZnYLsAKIBRY75zab2Y1N9y90zn1sZm8AG4AA8LhzblNbBheJNu/nlfDAyu1cOmkA38ocdOxfEDmCOXfkdHj7yMzMdNnZ2Z48tki4Kaqo5YLfvUdS13heueV0usbrrBzSMjPLcc5ltnSfXjUiHmvwB1jw/Dqq6vw8f0OGylyOm145Ih77zd+2s3rXIX5zxQSG90n0Oo5EMJ3LRcRDb28r4qG3dzIraxCXTBrodRyJcCp0EY/sK6vhjmW5nNI3kR9epItVyIlToYt4oL5p3tzXENDFKiRkNIcu4oFfrdhGzu5Sfj9rEkNTu3sdR6KERugi7Wzlxwd4dFU+V01L46IJ/b2OI1FEhS7SjgpKq7njhfWM6d+D71842us4EmVU6CLtxNcQ4Jbn1uEPOB6arXlzCT3NoYu0k/ve2Eru3jIevjKDwSndvI4jUUgjdJF2sGLzp/zhH7uYe2o6F4zr53UciVIqdJE2tudgNXe+uJ7xA3tyz4WjvI4jUUyFLtKGauv9fOe5HAx4aHYGCZ00by5tR3PoIm3oZ69tYVNhBY/NyWRQclev40iU0whdpI28klvIMx/uYf4ZQzl3dB+v40gHoEIXaQN5RYf57+UbyUxP4rvnjfQ6jnQQKnSREKvx+bn52bV0jovl97MnERert5m0D82hi4TYD17ZxPaiSp6al0W/nl28jiMdiIYOIiH0QvZe/phTwIKzhnHGiFSv40gHo0IXCZHN+8r5wZ82cdrJJ3HrOSO8jiMdkApdJATKa+q56Zm19Ooax+9mTSI2xryOJB2Q5tBFTpBzjjtfXM++shqWfXsaKd0TvI4kHZRG6CIn6NFV+by15QD3XDCKyenJXseRDkyFLnICPth5kPvf2MqF4/sx7/TBXseRDk6FLnKciipqWfD8OgandOO+y8Zjpnlz8Zbm0EWOQ70/wM3PraWqroHnbphK9wS9lcR7ehWKHIf739jKmk9K+e3MiYzok+h1HBFAUy4irfbn9ft47L1dzDk1nRkTB3gdR+RzKnSRVtj6aQXf++MGMtOTdJFnCTsqdJEglVfX8+2nc0js3ImHr8wgvpPePhJeNIcuEoRAwHHbsnXsK6th6fxp9O7R2etIIl8Q1BDDzKab2TYzyzOzu79kvSlm5jezy0MXUcR7D6zcwdvbirn3ojE6eEjC1jEL3cxigYeA84HRwCwz+8LkYdN69wErQh1SxEtvbTnA71bu4PLJA7lqaprXcUSOKpgRehaQ55zLd875gKXAjBbWWwC8BBSFMJ+Ip3YWH+aOZbmMG9CTn108VgcPSVgLptAHAHubLRc03fY5MxsAXAIs/LI/ZGbzzSzbzLKLi4tbm1WkXVXWNm4EjesUw8KrJ9M5LtbrSCJfKphCb2lI4o5YfgC4yznn/7I/5Jxb5JzLdM5lpqbq5P8SvvwBx21Lc9lVUsWDsyYxoJeuPCThL5i9XAqAQc2WBwL7jlgnE1ja9HE0BbjAzBqcc38KSUqRdvbrN7excmsRP5kxhtOGpXgdRyQowRT6GmC4mQ0BCoGZwOzmKzjnhnz2s5k9CfxFZS6R6pXcQh5+ZyezstK4elq613FEgnbMQnfONZjZLTTuvRILLHbObTazG5vu/9J5c5FIsn5vGd/74wayBifz42+M0UZQiShBHVjknHsdeP2I21oscufcNSceS6T9FVXUMv/pbFK6J/DIVToSVCKPjhQVAWrr/cx/OofK2gZeuuk0TtJl5CQCqdClw3POcc/yjeTuLWPhVRmM6tfD60gix0WfKaXDe+TdnSxfV8jt54xg+th+XscROW4qdOnQ/rJhH/e/sY1vTOjPgq8N8zqOyAlRoUuHlbO7lDteWE9mehL3Xz6emBjt0SKRTYUuHdLug1XcsCSb/j07s2hOpg7rl6igQpcOp6zax7wn1xBwjifmZZHcLd7rSCIhoUKXDsXXEODbT+dQcKiGRVdnMiSlm9eRREJGuy1Kh+Gc4+7lG1i96xAPXDGRrCG6UIVEF43QpcP4zd92sHxt4+6JF08acOxfEIkwKnTpEJ7+cDe/W7mDb04eyH+drd0TJTqp0CXqvb5xP/e+somzT+nNzy8dpxNuSdRSoUtU+2DnQW5bmktGWhIPzs6gU6xe8hK99OqWqLV5Xznzl2QzOKUrf5ibSZd47Wsu0U2FLlFpz8Fq5i5eQ2LnTjx1bRa9umpfc4l+KnSJOsWVdVy9eDUNgQBLrsuiX09dD1Q6BhW6RJXy6nrmLv6IAxW1LL5mCsN6J3odSaTdqNAlalTW1jPniY/IKzrMo1dnkpGW5HUkkXalQpeoUFXXwLwn1rC5sJyHr8zgqyNSvY4k0u5U6BLxauv9XP9UNmv3lPLbmZM4Z3QfryOJeELncpGIVtfg59tP5/DhroP837cmcOF4XXFIOi6N0CVi1fsD3PLcOt7dXszPLxnHJZMGeh1JxFMqdIlI9f4Aty3N5a0tB/jxN8YwMyvN60gintOUi0ScugY/C55bx5tbDvA/F4xi7mmDvY4kEhZU6BJRauv93PhMDu9sK+ZHF43mmtOHeB1JJGyo0CViVPsauP6pbD7IP8jPLx3HLE2ziPwbFbpEhMraeuY9sYa1e0r59TcncGmGNoCKHEmFLmGvvLrxCNDNheX8flaGdk0UOQoVuoS14so65i5uPJz/kasmc64OGhI5KhW6hK384sPMfeIjSip9PDY3U4fzixxDUPuhm9l0M9tmZnlmdncL919pZhuavv5pZhNCH1U6knV7Srl84QdU1fl5fv40lblIEI45QjezWOAh4FygAFhjZq8657Y0W20X8FXnXKmZnQ8sAqa2RWCJfis/PsDNz62ld2Jnnro2iyEp3byOJBIRghmhZwF5zrl855wPWArMaL6Cc+6fzrnSpsUPAe2CIMdl6Ud7uGFJNsN7J/LSTaepzEVaIZg59AHA3mbLBXz56Ps64K8t3WFm84H5AGlp2odY/sU5x29X7uCBv+3gjBGpPHJlBt0StIlHpDWCecdYC7e5Flc0O4vGQv+Plu53zi2icTqGzMzMFv+GdDy19X7+e/lGXl5XyGUZA/nFZeOIi9VphkRaK5hCLwAGNVseCOw7ciUzGw88DpzvnDsYmngS7Q5U1DL/6RzW7y3jjnNHsOBrwzBraQwhIscSTKGvAYab2RCgEJgJzG6+gpmlAcuBq51z20OeUqJS7t4y5i/J5nBdAwuvmsz0sX29jiQS0Y5Z6M65BjO7BVgBxAKLnXObzezGpvsXAvcCJwEPN42uGpxzmW0XWyLd8rUF3L18I70TE1h+3Wmc0reH15FEIp45581UdmZmpsvOzvbkscU7/oDjvje2smhVPtOGJvPwlZNJ7hbvdSyRiGFmOUcbMGs3Amk3RZW13L4sl/fzDnL1tHTuvWi0Nn6KhJAKXdrF+3kl3Lo0l8raeu67bBxXTNFuqyKhpkKXNuUPNO5f/vu/72BoSjeevX4qI/smeh1LJCqp0KXNHKio5dal6/gw/xCXZQzkpxePoWu8XnIibUXvLmkT72wr4v+9sJ5qn59ffXMCl0/W2SBE2poKXULqcF0D//vaxzz/0R5G9OnO0tkZDO+jKRaR9qBCl5D5584SvvfHDRSW1fDtM4Zy+7kj6BwX63UskQ5DhS4nrNrXwP1vbOPJf37CkJRu/PHGU5mcnux1LJEOR4UuJ2TNJ4f47ovr+eRgNdecNpi7pp9Cl3iNykW8oEKX43Koysd9f93Ksuy9DEruwtL505g29CSvY4l0aCp0aZVAwLF0zV7uX7GVw7UNzD9jKLeePVznLhcJA3oXStA2FpTz/Vc2sX5vGVlDkvnZxWMZoT1YRMKGCl2O6VCVj9+8tZ1nVu/mpG4J/OaKCVw8cYDOWy4SZlToclTVvgYW/2MXC9/Np9rXwJxp6dzx9ZH07BLndTQRaYEKXb6gwR9gWfZeHvjbDoor6zh3dB/umj6SYb01vSISzlTo8jnnHCs2f8r9K7aRX1zF5PQkHrkyg8zB2qdcJBKo0AV/wPHXTft56O2dfLy/gmG9u/PYnEzOGdVb8+QiEUSF3oHV+wO8vK6Qhe/sJL+kiqGp3fjVNydw8cT+dNKFJ0Qijgq9A6rx+Xkhey+LVuVTWFbD6H49ePjKDM4b05fYGI3IRSKVCr0D2XOwmmdW72bZmr2U19STmZ7Ezy4Zy5kjUjW1IhIFVOhRLhBwvLujmKc/2M3b24qIMeO8MX2Ye+pgsoYkq8hFoogKPUodqKjlT+sKef6jPXxysJqU7gks+NpwZmel0bdnZ6/jiUgbUKFHkWpfA29uPsBLawt4P6+EgIPM9CTu+PpIpo/pS3wnbegUiWYq9AhX7w+wOv8QL68r5I1N+6ny+RnQqws3nzWMSyYNYGhqd68jikg7UaFHoBqfn1U7ilmx+VNWflxEeU09iQmd+M/x/bk0YwBTBicTo71VRDocFXqEKK6sY9X2Yt7c8invbi+mtj5Aj86dOGdUH74+pg9njuyty72JdHAq9DBV4/OzetdB3s8r4b0dJWz9tBKAvj06863MQZw3pi9ZQ5KJ0wFAItJEhR4mKmvryd1bRs7uUlbnHyJndyk+f4D42BgyByfxvekj+cqwVMb076HpFBFpkQrdA4GAY9fBKjYUlJH9SSk5u0vZdqAS58AMTunbg7mnpfMfw1PJGpysa3SKSFBU6G2stt7P9gOVbN5XwZZ9FWzeV87WTyup9vkB6J7QiUlpvZg+ti+T05OYOKgXiZ11vnERaT0VeggEAo6iyjrySw6TX1zFzuJ/fS8sq8G5xvUSEzoxqn8PvpU5iNH9ezBuQE9G9EnU+VNEJCSCKnQzmw78FogFHnfO/eKI+63p/guAauAa59zaEGf1RCDgKKupp7iyjuLKOg5U1FJYVkNhaQ0FZdUUltawr6wWnz/w+e90jY9laGo3MtKSuHzyQEb2SWRM/54MTOqi+W8RaTPHLHQziwUeAs4FCoA1Zvaqc25Ls9XOB4Y3fU0FHmn67innHP6Aw+cPUOPzU+3zU+VroKrOT7WvgWqfn8raBsqqfZTX1FNWXd/4vaae0iofxZV1lByuoyHgvvC3UxMTGJjUhbEDejJ9bD8GJHVhyEndOLl3N/r26KxzpIhIuwtmhJ4F5Dnn8gHMbCkwA2he6DOAJc45B3xoZr3MrJ9zbn+oA7+zrYifvfYxgYAj4BwBBwHncK7xQg31/gA+fwBfQ+N398UublGMQY8ucfTqEkfPrvEkd4vnlL6JpCYmkJqYQO/Ezp//3K9nZ+3zLSJhJ5hCHwDsbbZcwBdH3y2tMwD4t0I3s/nAfIC0tLTWZgUgsXMcI/skYgaxMUaMGWYQY0aMQVxsDHGxMSR0iiG+U8zny13jY+kaH0u3hE7/9j0xIY6eXeNITOik6RARiWjBFHpLLXfkuDeYdXDOLQIWAWRmZgY5dv53k9OTmJyedDy/KiIS1YI5zLAAGNRseSCw7zjWERGRNhRMoa8BhpvZEDOLB2YCrx6xzqvAHGs0DShvi/lzERE5umNOuTjnGszsFmAFjbstLnbObTazG5vuXwi8TuMui3k07rY4r+0ii4hIS4LaD9059zqNpd38toXNfnbAzaGNJiIiraFT9YmIRAkVuohIlFChi4hECRW6iEiUMBfssfGhfmCzYmD3cf56ClASwjihEq65IHyzKVfrKFfrRGOudOdcakt3eFboJ8LMsp1zmV7nOFK45oLwzaZcraNcrdPRcmnKRUQkSqjQRUSiRKQW+iKvAxxFuOaC8M2mXK2jXK3ToXJF5By6iIh8UaSO0EVE5AgqdBGRKBHWhW5m081sm5nlmdndLdx/ipl9YGZ1ZnZnGOW60sw2NH3908wmhEmuGU2Zcs0s28z+IxxyNVtvipn5zezycMhlZmeaWXnT85VrZveGQ65m2XLNbLOZvRsOuczsu82eq01N/5bJYZCrp5n92czWNz1f7XI22CByJZnZy03vyY/MbOwJP6hzLiy/aDxV705gKBAPrAdGH7FOb2AK8L/AnWGU6zQgqenn84HVYZKrO//abjIe2BoOuZqt93caz+p5eTjkAs4E/tIer6tW5upF4zV905qWe4dDriPWvwj4ezjkAu4B7mv6ORU4BMSHQa5fAj9s+vkUYOWJPm44j9A/vzi1c84HfHZx6s8554qcc2uA+jDL9U/nXGnT4oc0XsEpHHIddk2vHqAbLVwm0ItcTRYALwFF7ZCpNbnaWzC5ZgPLnXN7oPF9ECa5mpsFPB8muRyQaGZG46DmENAQBrlGAysBnHNbgcFm1udEHjScC/1oF572WmtzXQf8tU0TNQoql5ldYmZbgdeAa8Mhl5kNAC4BFtJ+gv13PLXpo/pfzWxMmOQaASSZ2TtmlmNmc8IkFwBm1hWYTuP/oMMh14PAKBovi7kRuNU5FwiDXOuBSwHMLAtI5wQHf+Fc6EFdeNoDQecys7NoLPS72jRR08O1cFtLF+p+2Tl3CnAx8NM2TxVcrgeAu5xz/nbI85lgcq2l8bwZEwXYL44AAAH9SURBVIDfA39q81TB5eoETAYuBM4DfmBmI8Ig12cuAt53zh1qwzyfCSbXeUAu0B+YCDxoZj3CINcvaPwfcy6Nn1DXcYKfHIK6YpFHwvXC00HlMrPxwOPA+c65g+GS6zPOuVVmdrKZpTjn2vLkRcHkygSWNn4iJgW4wMwanHNtWaDHzOWcq2j28+tm9nCYPF8FQIlzrgqoMrNVwARgu8e5PjOT9plugeByzQN+0TTdmGdmu2ics/7Iy1xNr695AE3TQbuavo5fW2+0OIGNCp2AfGAI/9qoMOYo6/6I9tsoesxcQBqN11c9LZyeL2AY/9oomgEUfrYcDv+OTes/SftsFA3m+erb7PnKAvaEw/NF4/TByqZ1uwKbgLFe52paryeNc9Td2vrfsBXP1yPAj5p+7tP0uk8Jg1y9aNo4C9wALDnRxw3bEboL4uLUZtYXyAZ6AAEzu43GLckVR/3D7ZALuBc4CXi4adTZ4Nr4jG9B5roMmGNm9UANcIVrejV5nKvdBZnrcuAmM2ug8fmaGQ7Pl3PuYzN7A9gABIDHnXObvM7VtOolwJuu8dNDmwsy10+BJ81sI41TIXe5tv2UFWyuUcASM/PTuNfSdSf6uDr0X0QkSoTzRlEREWkFFbqISJRQoYuIRAkVuohIlFChi4hECRW6iEiUUKGLiESJ/w8hvp1TJ/0RGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "krw_plot = [krw(sw) for sw in sw_plot]\n",
    "kro_plot = [kro(sw) for sw in sw_plot]\n",
    "fw_plot = [fw(sw) for sw in sw_plot]\n",
    "plt.figure(1)\n",
    "plt.plot(sw_plot, krw_plot, sw_plot, kro_plot)\n",
    "plt.show()\n",
    "\n",
    "plt.figure(2)\n",
    "plt.plot(sw_plot, fw_plot)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the grid\n",
    "mesh = Grid1D(dx = L/nx, nx = nx)\n",
    "x = mesh.cellCenters\n",
    "\n",
    "# create the cell variables and boundary conditions\n",
    "sw = CellVariable(mesh=mesh, name=\"saturation\", hasOld=True, value = 0.)\n",
    "# sw.setValue(1,where = x<=dx)\n",
    "sw.constrain(1-sor,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some tests on fipy numerix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "krw_cell = krw(sw)\n",
    "krw_cell.value\n",
    "# It works fine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,\n",
       "       0.0046875, 0.0046875])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "eq = TransientTerm(coeff=1) + UpwindConvectionTerm(coeff = u\n",
    "*(sw**2./muw)/(sw**2./muw+(1-sw)**2./muo) * [[1]]) == 0\n",
    "\n",
    "sw.constrain(1.,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n",
    "\n",
    "steps = 100\n",
    "viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)\n",
    "for step in range(steps):\n",
    "    sw.updateOld()\n",
    "    swres = 1.0e6\n",
    "    while swres > 1e-5:\n",
    "        swres = eq.sweep(dt = dt, var = sw)\n",
    "        print(swres)\n",
    "        \n",
    "viewer.plot()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
